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Chapter | 
INTRODUCTION 
1.1 The Marine Riser Problem 
1.1.1 Description of the Marine Riser System 


A marine riser 1S a long, slender pipe used in offshore drilling 
operations. Extending from a fixed or floating platform at the sea sur- 
face to a wellhead connection at the seafloor, the riser contains the drill 
string, which it guides and supports, and the drilling mud, which it returns 
to the surface for reuse. Figure 1.1 depicts a typical riser installation. 

Risers with outside diameters of less than two feet are commonly 
used in water several hundred feet deep.!~3 Such a long, slender structure 
tends to buckle under its own weight unless it is supported in some manner. 
To prevent buckling and reduce deflections, a tension force, which is 
usually somewhat greater than the total submerged weight of the riser sys- 
tem, 1s applied at the surface end of the riser. 

Near the upper end of the riser, there is normally a slip joint 
which permits lengthening and shortening of the riser to compensate for 
changes in water depth caused by tides and waves. At the bottom end, where 
the riser connects to the blowout preventer, there is a flexible ball joint 
which permits some angular rotation without excessive bending when misalign- 
ment occurs. 

In addition to its own weight and the supporting tension, the 
riser is also subjected to the hydrodynamic forces of currents and water 
waves and to stresses resulting from the relative lateral displacement of 


the riser ends, which occurs when the surface platform drifts off station. 





1.1.2 The Riser Differential Equation 


Figure i.2 is a free body diagram of the deflected riser showing 
the coordinate system employed in this thesis and the loads involved in 
the problem. All applied loads and deflections are assumed to be coplanar. 
The governing differential equation, which is derived by applying 
the equations of static equilibrium to the element of riser length shown 


in Fig. 1,3, fs 


a, (EI fy - & (1(x)§4) = p(x) en 
where 

—E = the modulus ot elasticity of the riser material 

i = the moment of inertia of the riser cross section 


Tixj= the axia! tension in the riser 
p(xj- the lateral force per unit length of the riser 


Summaticn of forces in the vertical direction yields 


dix) : 
— f 
dx wi x) (lik) 
where 
w(x) = the submerged weight of the riser system per unit length. 


For risers considered herein, E, I, and w are constant and Eq. 1.1 reduces 


to 
ae dé d 
EL - Tix) =$- woe = p(x) (153) 
dx dx 5 


When the lateral force tntensity is time dependent, the response is also 


time dependent, and the equation of equilibrium 1s 


2 4 2 
ory y oy Ls a OY . ae (1.4) 
m ae: COO) att EI m We ae We pO) 


in which m is the riser system mass per unit length and c(x) is a damping 
function. 

in a structural sense, the marine riser is simply a tension beam 
whose axial tension varies along its length. Because of the varying ten- 
sion, the governing differential equation does not possess a closed form 
solution. Therefore, in the structural analysis of the marine riser, it is 


necessary to employ an approximate method of solution. 
1.2 Importance of the Problem 


In riser installations, it 1s necessary to limit the bottom rota- 
tion to a few degrees.*+* The maximum rotation must be less than that per- 
mitted by the design of the ball joint, which is usually about ten degrees. 
It is seldom possible, however, to attain such a rotation without producing 
excessive stresses in the riser. Moreover, in most installations, if the 
rotation is more than a few degrees, the drill pipe rubs against the riser, 
the ball joint, the blowout preventer, and the well casing resulting in 
wear which is great enough to constitute fa‘ lure. ° 

Bottom rotation 1s controlled, in part, by the positioning system 
used to keep the surface support platform on station. Whether a conventional 
mooring system or a dynamic positioning system is used, accurate information 
describing the behavior of the riser in its operating environment is needed 
for proper design of the system. 

Offshore drilling operations are expensive. The failure of a 
marine riser system results not only in the loss of costly hardware, but 


also in the loss of considerable productive time which is valued in thousands 


of dollars per day. Thus, there is much practical significance in any 
method of analysis which leads to a better understanding of riser behavior 
under random environmental loads and reduces the risk of unexpected failure. 

Published methods of riser analysis consider the effects of cur- 
rent, top offset, and top tension, but often ignore dynamic wave loads. 
Because the marine riser is 4 long, flexible structure, its fundamental 
period of vibration is typically on the order of a few seconds. Ocean 
waves genera'ly have periods ranging from six to thirty seconds. Hence, it 
is not uncommon for one or more of the riser natural periods to be in the 
range of wave periods, where the dynamic response is significantly greater 
than the static response. Neglecting the dynamic nature of the wave forces 
may resuit in unacceptable error in response calculations. 

Even when wave torces are considered, the approach used in pub- 
lished analyses is to calculate deterministically the response to a design 
wave Of a given height and frequency. Anyone who has seen the ocean will 
agree that unlike a simple design wave, real waves are extremely complex 
and highly irregular. A more realistic alternative to the design wave 
approach “s to treat the wave forces as the random physical phenomenon which 
they are and calculate the random response caused by these forces. Such 


an approach ‘s used in this thesis. 
1.3 Summary of Previous Work 


There are two separate areas of previous research which form a 
basis for this thesis. The first of these is the development of approximate 
methods of analyzing beams with varying axial tension. Fischer and Ludwig’ 
utilized the infinite power series method to calculate the response of a 


marine riser to a static current force and constant top deflection. This 


approach was also employed by Huang and Dareing°? to determine the natural 
frequencies of beams with variable tension. Frohrib and Plunkett® used 
a perturbation solution to determine the natural frequencies of a drill 
string. Finite difference methods were employed in the analysis of the 
driiling riser which was to have been used in Project Mohole’»® and have 
been uti!ized in the study of marine risers up to 1,300 feet long for 
offshore oi) drilling ‘2% Because the finite difference method permits a 
simple transformation of the governing differential equation into matrix 
form, @ torm which is conventent for machine computation, that method is 
used here for the calculation of riser response to random environmental 
forces. 

The second body of previous research to which this thesis is 
indebted s the development of models which accurately describe the random 
nature of ocean waves and the random forces they cause on submerged objects. 
Excellent comprehensive summaries of this development have been given by 
Pierson-° and Neumann and Pierson.!? Some of the more noteworthy mile- 
stones will be mentioned here. Fundamental to this development was the 
work of Rice 2 on the statistical properties of Gaussian processes. The 
demonstration by Pierson and others'® that the sea surface is essentially 
a Gaussian random process paved the way for the applications of Rice's 
work to the study of ocean waves. Longuet-Higgins * showed that the random 
wave heights for a narrow band spectrum have a Rayleigh probability distri- 
bution. From this early theoretical work and oceanographic observations, 
several models were formulated for the spectral density function of the 
sea surface elevation. /4-!6 

A significant extension of the work on ocean wave spectra was 


accomplished by Borgman,!7»!8 who derived expressions for the spectral 


density functions of ocean wave forces on piles and demonstrated their 
validity by comparison with observed force spectra. BOrgman's work opened 
the door to treatment of the response of ocean structures to wave forces 

as a random vibration problem. Such an approach was used by Foster!? 

and Malhotra and Penzien,2°>2! who adapted Borgman's model to the nondeter- 
ministic analysis of offshore towers; such an approach is also used in this 


thesis in the nondeterministic analysis of marine risers. 
1.4 Objective and Scope 


The objective of this thesis is to develop a mathematical model 
for predicting the random dynamic response of a marine riser to a combina- 
tion of random wave forces, deterministic steady current forces, and random 
operational movement of the surface support platform. Wave induced platform 
motion is not included in the model. 

The finite difference method is used to transform Eq. 1.3 into 


the matrix form 


CKaxn (YF axl {Ph yyy (1.5) 
where 
[K] = a stiffness matrix whose elements are derived by applica- 
tion of difference equations at each of n nodes along the 
riser axis 
iy} = a vector of node deflections 
{pi = a vector of force intensities at the nodes 


While it is expected that the accuracy of the finite difference model depends 
on the number of nodes used in the model, there is no way to determine a 


priori how many nodes must be used to accurately represent the real structure. 





In Chapter 2, the effect of the number of nodes on model accuracy is investi- 
gated, and criteria are established for selecting the number of nodes. 
Treatment of the riser problem in this thesis differs from pub- 
lished methods of analysis in that the wave forces and top offset are con- 
Sidered to be Gausstan, stationary random variables. It is shown that the 
riser response is also a Gaussian, stationary random variable, whose proba- 
bility density function is thus completely determined by the mean and 
variance. Because the mean response can be set to zero by a transformation 
of coordinates, knowledge of the response variance is sufficient to make 
probability calculations. The response variance is obtained by integrating 
the response spectral density function with respect to frequency. Thus, 
One objective of this thesis is to develop expressions for the spectral 
density functions of response parameters which resuit from a random sea 
represented by the oceanographer's sea surface elevation spectrum. 
Initially, only the response caused by random waves is considered. 
In Chapter 3, static response spectral density functions are derived, and 
examples are given to show how nondynamic effects influence riser response. 
Dynamic response spectra! density functions resulting from waves alone are 
derived in Chapter 4. Modification of the hydrodynamic forces by wave 
structure interaction is considered. The model is completed in Chapter 5 
where dynamic response spectral density functions caused by a combination 
of random waves and a steady deterministic current are derived, and the 
influence of random top offset on response variance is determined. In 
Chapter 6, the results of several sample calculations are used to show the 


effects of various problem parameters on the response. 





1.5 Notation 


Each symbol used in this thesis is defined where it first appears 


in the text. For ease of reference, the notation is summarized here. 


[B] transformation matrix relating random moments to 
normal coordinates 


aS coefficient equal to Cy 0 ro" 

C. coefficient equal to 5 Cp o D 

Ch drag coefficient 

Cy inertia coefficient 

[C,] modal damping matrix 

Fe] diagonalized modal damping matrix 
D riser outside diameter 

E modulus of elasticity 

ell expectation operator 

2) energy density per unit sea surface area 
NE} error vector 

F probability distribution function 
Gy top tension ratio 

G. length ratio 


H(2) 


[J] 


[kK] 


Lim 


wave height 


complex frequency response function 


moment of inertia of riser cross section 


transformation matrix relating random moments to 
deflections 


riser stiffness matrix 


riser length 


4 


(n+1) 


parameter equal to 
sf 


wave length 


random bending moment in riser 


random hydrodynamic force intensity 


force intensity equivalent to top offset 


vector of random modal forces 


probabi'ity 


correlation function 


one-sided spectral density function 


two-sided spectral density function 


axial tension in riser 


tension at top end of riser 


U(x,t) 


U(x,t) 


V(x,t) 


W(x,t) 


News t) 


[c] 
[c], {c! 


imc ot 


tension at bottom end of riser 


TL? 


parameter equal to ae 
(nt+1)°EI 


random horizontal water particle displacement caused 
by waves 


horizontal component of random water particle velocity 
due to waves 


horizontal component of wave induced water particle 
displacement with respect to structure 


random vertical component of water particle velocity 
due to waves 


random horizontal deflection of riser due to waves 
and current 


random horizontal detlection of riser due to top 
offset 


total rendom horizontal deflection of riser 
random top offset 


vector of influence coefficients relating top offset 
to riser deflections 


structure damping matrix 
equivalent linear damping matrices 
environmental damping matrix 


base of natural logarithms 


en JC 


nc 

fees nk, nz 
nm 

p(x,t) 


ip 


probability density function 
gravitational constant, 32.2 ft/sec” 
water depth 

impulse response function 

/-1 

subscripts of spectral components 
wave number 

mass ot riser system per unit length 

etfective mass of riser system and surrounding water 
number of nodes in finite difference model 

number of spectral components 

subscripts on riser model nodes 

number of natural modes in solution 


laterai force intensity on riser 


vector of force intensities due to lateral loads and top 
offset 


spatial wave force modifier 


subscripts on natural modes 


time 





A 2 


duration of environmental conditions 

horizontal component of water particle velocity 
steady current velocity 

horizontal component of water particle acceleration 
wind velocity 


submerged weight of riser system per unit length 


wl? 


a parameter equal to ora: a 
2(n+1)~EI 


vertical coordinate 
horizontal coordinate 
power transfer function 
random bottom rotation 


ratio of finite difference solution to analytical 
solution 


summation operator 
velocity potential 
random phase angle 
spectral frequency 


ratio of current velocity to standard deviation of rela- 
tive water velocity with respect to structure 





Ig 


B internal friction damping ratio 

Y unit weight of water 

6 Dirac delta function 

C total damping ratio 

n(t) sea surface elevation 

u mean value; tension beam parameter equal to YT/EI 

v radius of gyration of area under spectral density 
function 

E dummy time variable 

) mass density of water 

Oo standard deviation 

a variance 

1 time lag 

> mode shape 

Ww natural frequency of riser 


Each dot over a symbol indicates a derivative with respect to 


time. Thus, U(x,t) is the horizontal component of random water accelera- 


tion caused by waves. 
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Chapter 2 
A MATHEMATICAL MODEL OF THE MARINE RISER 


2.1 General 


In this chapter, the governing differential equation of the marine 
riser is transformed into matrix form using the method of finite differ- 
ences. This mathematical model of the structure is then tested by compar- 
ing the results of finite difference analysis to known exact solutions and 
published approximate solutions, and the number of nodes necessary to accu- 
rately model the structure is determined. In the process of establishing 
the validity of the fin te difference structural model, a modified wave 


force mode! ©s aiso evolved and tested 
2.2 The Finite Difference Structural Model 


In the finite difference method, a continuous structure having 
an infinite number of degrees of freedom is represented by a mathematical 
model having a finite number of degrees of freedom. The differential equa- 
tion of the continuous structure is transformed into a system of simultan- 
eous algebraic equations by substituting difference equations for the de- 
rivatives in the differential equation at each of a finite number of discrete 
points or nodes along the riser length. 

Consider the riser to be divided into n+i increments of equal 
length, Ax = at joined by n interior nodes as shown in Fig. 2.1. Nodes 
0 and n+l are at the bottom and top ends of the riser respectively, and 
nodes -1 and n+t2 are imaginary nodes located a distance Ax below the bottom 


and above the top. The finite difference model is formed by writing 





Eq. 1.3 at each of the interior nodes and replacing the derivatives with 


standard central difference equations. 22 
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Multiplication of each term of Eq. 1.3 by 
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and Xn is the x coordinate of node ni. The terms of Eq. 2.6 may be con- 
veniently rearranged so as to combine coefficients of the node deflections. 
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When Eq. 2 8 is applied at each of the nodes, there result n simultaneous 
linear algebraic equations in n+4 unknown deflections. 

The number of unknown deflections is reduced to n by application 
of the four boundary conditions. At the lower end of the riser, x = 0, the 


deflection and moment are zero. 








1 ia (2.9a) 

| dy, 

M(O) = - El 5 

dx 
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Therefore, 

yo = <¥ (2.9b) 


The top deflection, yi. ° Ytop? which is the lateral displacement of the 


surface support platform, 1s considered as an input parameter and therefore 
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not an unknown, and the top moment is taken as zero. 
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Therefore, the deflection of the imaginary node n+2 is 


eye eh ton on (2.10) 


Equations 2.8, 2.9, and 2.10 are conveniently expressed in matrix 


form, 
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The stiffness motrix, [K], is a f1ive-term banded matrix whose nonzero ele- 


ments are 
‘ * 
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Elements of the vector iy} are the node displacements, and {p} 1s a vector 


of load intensities in which 
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The riser axial tension, which is used to compute les 1s determined 


by integrating Eq. 1.2 The tension at coordinate x is 

ix) = Ty + wx 2.14) 
where T° the axial tension at the lower end of the riser, 1s given by 
tee 1, = wh (2eliey 
and T is the axial tension at the surface. 


Equations 2.11, 2.12, and 2.13 describe the finite difference model 
of the marine riser in matrix form. It is of interest to compare this model 
with the finite difference model used by Butler, et al.,’ for engineering 
studies of the Project Mohole drilling riser. The formulation of the Mohole 
model was somewhat different from the formulation of the model presented 


herein. For the Mohole riser, axial tension was considered to be constant 


between nodes, and changes in tension were concentrated at the nodes. The 
Mohole mode! was also derived for the general case of variable EI, which 
could have been done here by applying the finite difference equations to 

Eq. I|.1 instead of Eq. 1.3. When the Mohole riser equations are written for 


the case of constant El, they are identical to those derived herein. 
2.3 Validity of the Finite Difference Structural Model 


In subsequent chapters of this thesis, the finite difference struc- 
tural model is used as a component of a larger mathematical model for calcu- 
lating the riser response to random wave forces. The validity of this larger 
mathematica! model depends, in part, on the validity of the finite difference 
structura! model. Therefore, before the finite difference model is applied, 
¥t 18 prudent to consider the question, "How well does the model represent 
the continuous ¢iser system?" 

In judging the validity of the finite difference model, two cri- 
teria are of particular interest. Because normal mode superposition is used 
in the dynamic analysis of the riser, its natural frequencies and mode shapes 
must be determined. The first criterion, then, is that the natural frequen- 
cies and mode shapes of the finite difference model must be reasonably close 
to the actual naturai trequencies and mode shapes of the continuous riser 
structure. To be effective, the finite difference mode! must also give node 
deflections, rotations, and moments which do not differ significantly from 
the actual deflections, rotations, and moments. Therefore, the second cri- 
terion 75 that the deflected shape of the finite difference model must reason- 
ably approximate the deflected shape of the actual riser. Since the random 
forces in the riser problem are those produced by ocean waves, it 1s espec- 


ially important to study the accuracy of the deflected shape of the model 
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when the applied force distribution is that which is typical of ocean 
waves. 

{9 order to form a judgment as to how well the finite difference 
model meets these two criteria, a standard for comparison is needed. The 
obvious standard 1s an analytical solution, 1f it exists. However, the 
marine riser has no closed-form solution, and therefore some other standard 
must be sought- 

One approach to testing the validity ot the finite difference 
model ts to perform the tests on a structure similar to the riser, but 
having an analytica! solution which can be used as the standard for compari- 
son. This approach is used in Section 2.4, where a constant tension beam 
1s modeled by the method of finite differences. The validity of this model 
is tested by compaing the natural frequencies, mode shapes, and deflected 
Shapes of the model with the corresponding analytical solutions to the 
constant tension beam. in the process of testing the deflected shape of 
the mode} under wave ‘oads, an improved wave force model is developed and 
tested. 

While validity tests of the finite difference method with the 
constant tension beam as a test structure are certainly worthwhile, some 
tests of the riser finite difference model are also definitely in order. 
in the absence of a closed-form solution to the riser equation, the only 
standard for comparison is the results of other approximate methods of 
analysis. In Section 2 5, the natural frequencies and deflected shapes of 


the riser fin te difference model are compared to other published solutions. 
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2.4 Tests of the Finite Difference Model of a Constant Tension Beam 
2.4.1 Formation of the Finite Difference Model 


The constant tension beam, which is shown in Fig. 2.2, is identical 
to the marine riser with one exception. Because dT/dx = 0, the governing 


differential equation of the constant tension beam is 


We 2 
ELS. EY = p(x) (2.16) 
dx dx 


it follows that the stiffness matrix for the finite difference model of the 
constant tension beam is identical to that given by Eqs. 2.12 for the riser 


except that Te 1s a constant and 


, ia 2 7 x. f 
Snisnt-1 > Kn niet a en 


_ 


2.4.2 Tests of Natural Frequencies and Mode Shapes 
Finite Difterence Solution 


The natural frequencies and mode shapes of the finite difference 
model are found by letting the force intensity vector of Eq. 2.11 represent 
the d’Aiembert forces, all other externally applied forces being zero. 

a ne 7 | 

{p} = -m pee (2-16) 


in which m is the mass per unit length of the riser system. Application of 


the method of separation of variables leads to the characteristic equation, 


[K] ty! = mw@ty) (2.19) 


solutions to Eq. 2.19 are the natural frequencies, TR ls 25 eee 


and mode shapes, roll}, of the finite difference model. 
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Analytical Solution 


The analytical solution for the natural frequency and shape of the 


h : F 
rt mode of a constant tension beam is23 


Van EI Te 
w >t = a= + 5 (2.20) 
r ic m ren em 
(r)yoy - rmx 
© (x) = sin ee (22210) 


It is convenient to use the following two dimensionless parameters in this 


study as well as tn the riser analysis which follows. Let 


oS 
et = ; 
G ce the length ratio 
and 
Ty o ° 
Gr = = the top tension ratio 


Then Eq. 2-20 may be rewritten 
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Equation 2.22 is the standard against which the natural frequencies of the 


finite difference model are to be compared. 
Comparison of Finite Difference and Analytical Solutions 


The natural frequencies and mode shapes of the finite difference 
model of the constant tension beam were calculated by using the University 


of {llinois Department of Computer Science System 360 Library routine EIGENP. 


(6, 


This routine computes the eigenvalues of a real, general matrix using the 
QR double-step method and computes the eigenvectors by inverse iteration.2* 

Comparison of the finite difference and analytical solutions was 
made for finite difference models having 3, 5, 7, 11, 15, 19, 23, and 31 
equally spaced interior nodes. In all cases, the first five modes were 
compared except for the model having only three interior nodes and thus only 
three natural modes. A top tension ratio of 0.5 was used to approximate a 
riser whose average tension is one-half of the total submerged weight. 
Length ratios of 8 and 1000, which are typical of risers located in shallow 
and moderately deep water, were used in the comparison. The accuracy of the 
finite difference model is expressed by the frequency ratio, which is simply 
the natural frequency of the finite difference model divided by the actual 
natural frequency. 

The results of the frequency comparison are shown in Fig. 2.3. As 
expected, the accuracy of the finite difference model improves as the number 
of nodes in the model increases. Also, the accuracy of the model is better 
for lower mode numbers than higher ones. It is significant that the finite 
difference model frequencies are more accurate for high G than for low G « 
As the length ratio increases, the structure becomes less like a beam and 
more like a vibrating string, which 1s apparently more suitable for finite 
difference analysis than a beam. Because G, 1S proportional to the water 
depth cubed, the finite difference model does a better job of predicting 
frequencies in deep water than it does in shallow water. Finally, although 
the actual number of nodes which one uses depends upon the accuracy desired, 
it is significant that the accuracy of the first three frequencies is 99 per- 


cent or better for the 31-node model. 
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The agreement between the mode shapes of the finite difference 
model and the analytical! solution is even better than the agreement between 
the frequencies. For each of the models tested, the first five mode shapes 
were normalized and the ratio of the deflection at each node to the corres- 
ponding analytical deflection was calculated. In all the models this ratio 


was 1.000 at every node. 
2.4.3 Tests of Deflected Shape Under Wave Loads 
Wave Force Distributions 


The second criterion for judging the validity of the finite dif- 
ference modei is that the deflected shape of the model under a given loading 
must reasonably approximate the actual deflected shape of the structure 
under the same loading. Because the model will be used to calculate the 
riser response to wave forces, a logical force distribution for the validity 
tests is one which is typical of water waves. The expression commonly used 


for the wave force on a vertical circular cylinder is Morrison's formula,?° 





p(x,t) = Cp 5 o Dulu| +C, p no" u (2.23) 
where 

p(x,t) = the wave force per unit length 

Cy = a dimensionless drag coefficient 

Cy = a dimensionless inertia coefficient 

p = the mass density of the water 

D = the outside diameter of the riser 

u(x,t) = the horizontal component of water particle velocity 


u(x,t) = the horizontal component of water particle acceleration 
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According to Eq. 2 23, the hydrodynamic force consists of two components, a 
drag force which is a function of water particle velocity and an inertia 
force which is a function of water particle acceleration. 

Water particle velocities and accelerations for use in Eq. 2.23 
are given by several wave theories found in the literature. Linear or Airy 
theory, which is derived for small amplitude waves, is used in this study. 
According to linear wave theory,¢© the horizontal component of water parti- 


cle velocity 1s 


Pel 7, SOK COSI KX) 7 
uo = ZH AG st Le sin (ky - at) (2.24) 


wn which 
H =~ the wave height measured from crest to trough 
g = gravitational acceteration (32.2 ft/sec‘) 
k = the wave number equal to 2a 
L = the wave length 
h = the water depth 
@ = the circular frequency of the wave 


in linear theov-y, convective acceleration 1s neglected as small, and the 
water particle acceleration is the time derivative of the water particle 


velocity 2° 


een *o07 =. 1 cosh (kx 7 
u — - > Hgk cosh h cos (ky at) (2.325) 


Of immediate interest is the spatial variation of the drag and 
inertia forces. From Eqs. 2.25 and 2.23, it is evident that the spatial 


variation of the inertia force is given by 


: _ cosh = 2.26 
a\x) cosh (kh ( 
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where q(x) is the ratio of the force intensity at coordinate x to the force 
intensity at the surface. Similarly, from Eqs. 2.24 and 2.23, it follows 
that the spatial variation of the drag force is 

q(x) = cosh (kx) (2.27) 

cosh” (kh) 

The spatial variation of inertia and drag forces for representa- 
tive values of kh is shown in Figs. 2.4 and 2.5. This spatial variation is 
completely defined by the parameter kh, which may assume a range of values. 
As the water depth decreases, kh approaches zero and the forces approach a 
uniform distribution. An upper limit on wave number may be determined with 
the help of linear wave theory which provides the following relationship 


between wave circular frequency, wave number, and water depth. 


9% = gk tanh (kh) (2.28) 


For large values of kh, the hyperbolic tangent approaches unity and 
2 
Q 
k Q (2.29) 
g 


With a wave circular frequency of 7/3 radians per second, which corresponds 
to a six second wave period, Eq. 2.29 yields a wave number of .034 feet” |. 
Therefore, an upper limit on kh is .934 h, where h is in feet. 
Comparison of Finite Difference and Analytical Solutions for Wave Inertia 
Force Deflections 

In this section, the finite difference deflections are compared 
to the analytical deflections of a simple beam subjected to the force given 
by Eq. 2.26. The problem is sketched in Fig. 2.6, and the analytical solu- 


tion is derived in Appendix A. Finite difference model deflections are 
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obtained by premu'tipiying both sides of Eq. 2 11 by cK. 


a= | ; 
‘yt = [K) {pi (2.30) 


in which 
cosh (kx, ,) 


Poi ~ Cosh (kh) (293) 
and [K] is the stiffness matrix of the constant tension beam with i. set to 
zero. A measure of the accuracy of the finite difference model is Aye the 
ratio of the node deflections computed by finite differences to the corres- 
ponding anaiytical deflections. 

The variation of midspan deflection ratio, M(L/2)? with kh is 
Shown tn Fig. 2.7 tor finite difference models having 3, 7, 15, and 31 equally 
spaced nodes. (it is helpful to remember that the number of length increments 
is one more than the number of nodes.) For values of kh less than 3, the 
finite difference model yields results which are in good agreement with the 
analytical solution. [In this region, errors are conservative and tend to de- 
crease as n increases. As kh increases above 3, the finite difference mode! 
becomes increasingly less accurate, although accuracy is improved somewhat by 
increasing n. Even with as many as 31 nodes, the error of the finite dif- 
ference mode’ ig greater than 10 percent tor kh greater than 37. In essence, 
the model acts as a tilter which removes those loads for which kh is large. 

If the model is to be useful in deep water where large values of 
kh occur, then it 1s obvious that a very large number of nodes must be used 
to accurately depict the force which, for large values of kh, decays rapidly 
with distance below the surface. For economy of computation, it is neces- 
sary to limit the number of nodes Thus, it is evident that an alternate 


force mode! is needed--one which accurately predicts the deflected shape, 
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even for large kh, without requiring an excessively large number of nodes. 


Such a model is developed in Section 2.4.4. 
2.4 4 The Lumped-Smeared Force Model 


A force mode’ which does not require a large number of nodes for sat- 


isfactory accuracy 1s the lumped-smeared force model shown schematically in 


Fig. 2.8. The lumped-smeared force intensity, Qni? is obtained by calculating 


1 
concentrated reactions to the distributed force, q(x), at each node and then 
distributing or smearing the total concentrated reaction at each node over the 
region halfway to adjacent nodes. 

Consider a section of the structure between nodes ni-] and ni+] shown 


in Fig 2.9a. Let Q and Q be the concentrated reactions at node ni. 


n14R ni,- 
Summatton of moments about node ni-1! in Fig. 2.9c yields 


ni 
el 7 of 5 Xni-1) a(x) dx (2.32) 


Xni-] 


Similarly, summation of moments about node nit! in Fig. 2.9d results in 


n+] Amit] 
nig . a (xps47 7X) a0) dx (2.33) 


Xni 


The lumped-smeared force intensity at node ni is, by definition, 


20 Xni+l 
= nv] + = bell 
Gni ape nae QnigRe h ] if [X44] 
x 


ni 
xn 
- x] q(x) dx +f [x - Xie q(x) dx (2.34) 
tie! 


It is evident that the lumped-smeared model preserves the total force on 


the structure. 
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Lumped-Smeared [Inertia Force 


The iumped-smeared inertia force intensity at node ni is found 
by substituting Eq. 2.26 into Eq. 2.34. Integration of the resulting ex- 


pression leads to 


2 
ee a [cosh (kx... 


= 2 COSh (Ke 
(kh a! \ Xni) 


g 
am cosh (kh) ni] 


+ cosh (kx - )j (2.35) 
ni*] 


Lumped-Smeared Drag Force 


To obtain the lumped-smeared drag force intensity at node ni, 
Eq. 2.27 1s substituted intc Eq 2.34. After integration and simplification, 


there results 


i ——___ ag Lcosh (2kx 


) 
8 cosh¢ (kh) ni-l 


- 2 cosh \2kx_- ) + cosh (2kx, 47) J + 4 (2.36) 


Tests of Deflected Shape Using Lumped-Smeared Wave Force Model 


The enelytical solutions for the deflections and end rotation of a 
constant tension beam under the inertia and drag forces of Eqs. 2.26 and 2.27 
are derived in Appendix B. The corresponding finite difference solutions 
are obtained from Eq. 2-30 with the stiffness matrix defined by Eqs. 2.12 
and Z 17 and elements of the force intensity vector given by Eqs. 2.35 and 
2.36. The validity ot the combined finite difference structural model and 
the lumped-smeared force model 1s measured by the parameter A, defined as 


the ratio of node deflection or rotation calculated with the model to the 


30 


corresponding detection or rotation given by the analytical solution. Hence, 
ey is the rotation ratio at x = 0 and A (L/2) is the midspan deflection ratio. 
In comparing the two solutions, the problem parameters have been 
chosen so that the constant tension beam approximates, as nearly as possible, 
typical! riser instailet'ons in the ocean. Accordingly, the load distribution 
parameter, kh, is varied from .01 to 10”, the former value corresponding to 
a practical’y unitorm toad and the jatter corresponding to a practically con- 
centrated load wery near the right end of the beam. Appropriate values of 
the parameter ., which appears in the analytical solution, are selected 
by noting that » may be written as a function of the riser parameters Gy and 
Gr. 


| pear emaraa 
ea il! 
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The length ratio, G > ot the constant tension beam is varied from 8 to 8000, 
é@ range which includes length ratios of typical riser installations. The top 
tension ratio, Gs, of a riser ts usually one or slightly larger. For this 
constant tension beam study, a ratio of r = 0.5 is used to approximate the 
average tension in a riser with a top tension ratio of 1.0. 

Some results of the compar?son between the finite difference solu- 
t?on with the lumped-smeared force model and the analytical solution are shown 
in Figs. 2.10 through 2.13. Figures 2.10 and 2.11 show the variation of M, 


0 


and A with number of nodes in the finite difference model for the 


We) 2) 
inertia force distribution of Eq. 2.26, for length ratios of 8, 216, 1000, 
and 8000, and tor load distributions parameters, kh, of 0.1, 10, 100, and 

100,000. The same information is shown in Figs. 2.12 and 2.13 for the drag 


force distribution of Eq 2.27. 
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As expected, the accuracy of the finite difference model with the 
lumped-smeared force system improves when the number of nodes is increased. 
Better accuracy with more nodes was also shown in Fig. 2.7 for the finite 
difference structural model without the lumped-smeared force model. Unlike 
Fig. 2.7, however, Figs 2.10 through 2.13 show no significant loss of ac- 
curacy as the parameter kh increases. Even when the applied force is essen- 
tially a concentrated force very close to the end of the beam, as it is when 
kh 1s TOP, the finite difference, lumped-smeared model yields results which 
are accurate to within one percent when ten or more nodes are used. Thus, 
the lumped-smeared force model corrects the main disadvantage of the unmodi- 
fied finite difference model--its inability to accurately represent the force 
distribution for al! values of kh. 

The lumped-smeared model has another significant characteristic. 
Except for small values of kh, as would be encountered in shallow water, 
the model displays a marked increase in accuracy as length ratio increases 
even when only a few nodes are used. This is particularly evident in the 
case of the end defiection ratios for G = 1000 and 8000. Even for a model 
with three nodes, when Gy = 8000 and kh is equal to or greater than 10, the 
error in end rotation is very small. Because the length ratio increases 
with the water depth cubed, the finite difference model with the lumped- 
smeared force system becomes increasingly accurate as water depth increases. 

With the assumption that the validity of the finite difference, 
lumped-smeared model for the constant tension beam is indicative of its 
validity for the marine riser, it may be concluded that the model is an 
acceptable method of representing the effect of wave drag and inertia forces 


on the riser structure. 
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2.5 Comparison of Riser Finite Difference Model Results With Other 
Approximate Solutions 
Although there is no closed-form riser solution to which the riser 
finite difference model results may be compared, there are published approxi- 
mate solutions which will serve as a standard for comparison. The literature 
includes approximate solutions for the natural frequencies as well as the 


static deflections of a marine riser. 
2.5.1 Comparison of Natural Frequencies 


The natural frequencies and mode shapes of the finite difference 
model of the marine riser were found by solving Eq. 2.19 with [K] given by 
Eqs. 2.12 Calculations were made by using the routine EIGENP which is 
described in Section 2.4.2. 

Huang and Dareing> have solved the characteristic differential 


equation of the marine riser, 


4 | 
dx ne 


using the infinite power series method, with y expanded into a power series 
in. X 


- Jc, x (2.38) 
jeo ° 
They calculated the first three natural frequencies of pipes having linear 
tension variation and several end conditions, one of which corresponds to a 
marine riser with zero bottom tension. Table 2.1 is a comparison of the 


first three eigenvalues as given by Huang and Dareing and as calculated by 


the finite difference model with thirty-one equally spaced nodes. In all 
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cases, the finite difference model frequencies are slightly lower than, but 
very close to, the series solution. 

Another approach to the solution of Eq. 2.37 is used by Frohrib 
and Plunkett® who employ a perturbation method. Again, a power series is 
used, but instead of powers of x, the series is in powers of the perturba- 


tion parameter, 


Re 
we 


Solution is obtained by setting coefficients of like powers of a equal to 
zero. 

Figure 2.14 is a comparison of the first three natural frequencies 
cf a marine riser with zero bottom tension as calculated by the finite dif- 
ference method, the series method of Huang and Dareing, and the perturba- 
tion method of Frohrib and Plunkett. Twenty-three equally spaced nodes 
were used for the finite difference solutions. The finite difference model 
and the series method give virtually identical results for all length ratios 
from zero to 1000. The perturbation solution agrees with the other two solu- 
tions for length ratios greater than 200. For lower values of Gy» the per- 
turbation parameter used by Frohrib and Plunkett increases rapidly, as 
Shown in Fig. 2.15. As in any perturbation solution, the results become 
less accurate as the magnitude of the perturbation parameter increases. 

The favorable agreement of the finite difference model frequencies 
with those of the series and perturbation solutions suggests that the former 
method is an acceptable one, at least for the first few modes of vibration 


and the range of length ratios studied. 
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2.5.2 Comparison of Responses to Steady Current Forces 


Fischer and Ludwig* solved Eq. 1.3 using the infinite power series 
method and presented their results in dimensionless form in a series of 
figures. From these figures may be calculated the top and bottom rotations, 
maximum moment, and top shear for a given top offset and two load distribu- 
tions. These two load distributions correspond to a uniform current and a 
linearly varying current with bottom current velocity equal to one half of 
the top current velocity. 

For the purposes of comparing the finite difference model with 
the series solution of Fischer and Ludwig, a test riser was selected having 


the following characteristics. 


D = 20 Inches 

Wall thickness = QO 438 inches 

L = 600 feet 

E = 29 x 10° pounds per square inch 
W = 248 pounds per foot 


Five combinations of top tension, lateral current load, and top 
oftset were investigated. Solutions for each of the five cases were obtained 
both by using the figures of Fischer and Ludwig and the finite difference 
model with various numbers of nodes. Two finite difference models were con- 
sidered, one with equally spaced nodes and the second with unequally spaced 
nodes. In the latter model, the node spacing in the center half of the riser 
is twice as long as the node spacing in the end quarters. The results for the 
five examples are given in Table 2.2. 

The two finite difference models yield results which are in close 


agreement with the series solution of Fischer and Ludwig. The agreement 


35 


improves as the number of nodes in the finite difference model increases, 


and even with as few as seventeen nodes, the agreement is very close. 
2.6 Selection of Number of Nodes in the Finite Difference Model 


The number ot nodes used in the finite difference model should 
strike a balance between economy and accuracy. The use of more nodes than 
are needed tor the accuracy desired wastes costly machine computation time, 
but the use of too few nodes may produce unacceptably inaccurate results. 

While this thesis does not profess to give detinitive criteria for 
selecting the number of nodes, the results of Sections 2.4 and 2.5 do provide 
useful guidance in tre absence of better criter°a. Because the lumped- 
smeared force moge! so effectively represents wave force distributions, the 
controlling criteria for selection of the number of nodes will usually be the 
accuracy with which the natural frequencies must be obtained. If the riser 
length ratio and the number of natural modes to be considered is known, 

Fig. 2.3 may be used to choose the number of nodes which correspond to the 
desired accuracy of the natural frequencies. Figures 2.10 through 2.13 
should also be checkea to insure that the accuracy of the lumped-smeared 


force model is also acceptable for the number of nodes selected. 
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Chapter 3 
THE STATIC RESPONSE OF A MARINE RISER TO RANDOM WAVE FORCES 


3.1 General 


In Chapter 2, a mathematical model of the marine riser was de- 
rived by applying the method of finite differences to the riser differen- 
tial equation, and a lumped-smeared mode! was developed to represent the 
wave force intensities at the nodes of the finite difference model. In 
this and subsequent chapters, these models are used to determine riser re- 
sponse to random wave forces. Consideration of dynamic effects is deferred 
until Chapter 4, so that in this chapter attention may be focused on the 
nondynamic effects which influence riser response. 

It wili be shown that the random wave force may be considered to 
be a stationary, zero mean, Gaussian random process. Riser response param- 
eters, the bottom rotation, deflections, and bending moments, which are 
linear functions of the wave force intensity, are also zero mean, Gaussian 
random processes because the class of Gaussian random variables is closed 
under linear operations.*° Thus, while it will not be possible to deter- 
mine the exact response of the riser at a particular time, the model should 
make it possible to calculate the probability that the response will] not 


exceed some tolerable ievel. 
3.2 Method of Solution 


As a typical response parameter, consider the random bottom rota- 


tion, Oy (t). (In this thesis, a random variable is represented by a capital 
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letters and the corresponding lower case letter denotes the state variable 
or range variable of the random variable). The probability that the bottom 


rotation does not exceed some allowable value, 6,° 1S given by 


P ure | oe j f. (6)de (3.1) 


In which t. 1s the probability density function. For a Gaussian bottom 
fe) 


rotation, the probability density function is 


(0 - u, )? 
fy (6) = exp - (3 2) 
re) Vv 21 o 25 
fe) | ) 
where 
Oo ae the variance of the bottom rotation 
0) 
O5 = the standard deviation of the bottom rotation 
fe) 
Uo = the mean value of the bottom rotation 
0 


From Eq. 3.2 it 18 seen that the probability density function of a Gaussian 
random variable is completely defined by its mean and variance. The mean 
values can be conveniently set to zero by a transformation of coordinates. 
Therefore, if it is possible to catculate the variance of the response, 
then it is also possible to determine the probability associated with a 


Specific level of response. 





Not all capitalized variables are random variables, however. For example, 
I, the moment of inertia of the riser cross section and E, the modulus of 
elasticity, are considered to be deterministic. 
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By definition, the variance is the second central moment of a 
random variable. 
o = Ele, - uy I (3.3) 
fe) fe) 
where E [] denotes the operation of taking the mathematical expection. 


For a zero mean random variable, the variance is the mean square value, 


ay |. 2 
oe, 7 & [95 ] (3.4) 


Suppose that a time record of the bottom rotation, Q,(t), were 
available. Because the variation of bottom rotation with time is a random 
process, this record would be just one realization of an infinite collection 
of possible records. This collection is called the ensemble, and the single 
record 1s calied a sample function. 

One of the statistical quantities which could be calculated from 
the ensemble 1s the autocorrelation function, Re which, by definition, 


is given by 


Roe (t, sty) = £E [9, (t,) e) 


» (tp) (3.5) 


The autocorrelation function is the ensemble average of the product of the 

random response at times t, and ty If the random process is stationary, 

then the autocorrelation function is inveriant with a shift in time. 
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Equation 3.6 states that the autocorrelation function is independent of 
time and depends solely on the time lag, 7. If the time lag is set to 
zero and the random process has a zero mean, then the autocorrelation func- 


tion 1s identically equal to the variance. 


"8 (0) = E [eo ,(t) a(t) ] = 2 (S57) 


Thus, 1f the autocorrelation function is known, the variance, and hence 
the probability density tunction, may be determined. 

The autocorrelation function is defined in the time domain. It 
is often convenient to shift from the time domain to the frequency domain 
by means of a Fourier transformation. By the Wiener-Khintchine theorem, 2’ 
the result of a Fourier transformation of the autocorrelation function of 


a stationary random process is the spectral density function, S, , (2), 
0-0 


x 


Le -i0T 4 (3.8) 
oa ] Ryo, (1) et ar 


where 2 1S the trequency in radians per second and i = v-1. Conversely, 
the autocorrelation function is the Fourier transform of the spectral 


density function. 


5 peolee Dare (3.9) 
0° “ 0°0 


Equations 3.7 and 3.9 suggest an alternate method of calculating 


the variance. 


2 720 
O = R (0) = f Ss. a ore dQ 
95 05% _ 99% 
ee (a) do (3.10) 
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The variance is seen to be the area under the spectral density function of 
a zero mean stationary random process. Thus, if the response spectral den- 
Sity function 1s Known, then the variance may be determined. 

Throughout this thesis, S(2) 1s used to denote a two-sided Spectrum, 
which 1s defined for both positive and negative frequencies and is an even 
function of frequency. An equivalent way of representing the same spectrum 
1s aS a one-sided spectrum, S(a) = 2 S(a). The one-sided spectrum, which 
1S nonzero for positive frequencies only, 1S more etficient for numerical 
1ntegration. 

In Sections 3.3, 3.4 and 3.5, a mathematical model 1s developed 
for determining the static response spectral density functions, $0,0,(%) 
ee and Su aM The subscript ni refers to node ni of the struc- 
tural model, Y(t) is the random deflection of node ni, and M(t) is the 
random moment at node iii The input to the model 1s the oceanographer's 
description of the random sea surface, the one-dimensional sea surface 
elevation spectral density function, Spt@)s where r{t) is the random el- 
evation of the sea surtace above mean water ieve!, as shown in Fig. 3.1. 
(Although ;, 1s a random variabie, 1t 1S represented by a lower case sym- 
bol in this thesis in order to agree with the notation commonly used in 
the literature). The static response spectra! density function is ob- 
tained from the sea surface elevation spectral density function through a 
series of intermediate relationships. These intermediate steps, which are 
developed next, are 

1. Calculation of the fiuid kinematic spectral densities from 


the sea surface elevation spectral density by means of linear 


wave theory, 


4] 


2. Calculation of wave force spectral densities from fluid 
kinematic spectral densities by means of the Morrison force 
formula, Eq. 2.23, and the lumped-smeared force model de- 
veloped in Chapter 2. 

3. Calculation of the static response spectral density from 
the wave force spectral density by means of the finite 
difference model developed in Chapter 2. 

Steps I and 2, which are treated in Sections 3.3 and 3.4, are based on the 
work of Borgman.:72!8 Because details of the derivation are not widely 


available, they are included here. 
3.3 Fluid Kinematic Spectra 


The four fluid kinematic spectral density functions which are 
used in the calculation of wave force spectral densities are 


Ss; e (2) 
Unk ne 


where 
Up, (t) = the horizontal component of the random water particle 
velocity at node nk 

U , it) = the horizontal component of random water particle ac- 


celeration at node nk 
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It is only necessary to derive an expression for S,, * (2) as a function 
nk~ nd 
of the sea surface elevation spectrum, Sin (@) because, as will be demon- 


strated, the remaining three spectra are readily obtained from Sy) U (aie 
nk~n2 


In this thesis, the sea surface elevation spectral density 
function is one-sided and defined in such a way that the area under the 
Spectrum is equal to the variance of the sea surface elevation. In the 
literature, there appear several formulas for the spectral density of either 
the wave height or the sea surface elevation for fully developed seas. Be- 
cause these spectra are not all defined in the same manner, one must use 
extreme care in using them for engineering analyses. Three of the more 


commonly used sea surface elevation spectra are described in Appendix C. 
3.3.1 Derivation of Expressions for Random Water Velocity and Acceleration 


For deterministic analyses, linear wave theory is often used to 
calculate the fluid velocities and accelerations associated with water 
waves. For the nondeterministic model developed here, linear wave theory 
is used to relate the fluid kinematic spectra to the surface elevation spec- 
trum, Sits In order to apply linear wave theory, it is first necessary 
to establish a relationship between the random sea surface elevation, n(t), 
and the sea surface elevation spectral density, San’ Such a relationship 
may be found by using energy considerations and the spectral decomposition 
theorem. 

The average total energy per unit surface area of the sea, Eo 
may be determined from the sea surface elevation spectrum, San’ It can be 


shown that2® 


=) of (3.11) 
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where 
y = the unit weight of the water 


o = the variance of the sea surface elevation 


Because the mean value of n is zero, by definition, the variance is the 


area under the spectral density curve. 


g = 
oO f Se (2) do (3,12) 
fe) 
Combining Eq. 3.11 and 3.12 yields 
ee ae s d 
: a nn (2) do (3.13) 
0 


Consider a typical sea surface elevation spectrum shown in Fig. 3.2. 
Let the spectrum be divided into nc frequency components, each of size Aa. 
It is assumed that the sea surface elevation is a stationary random process, 
which, by the spectral decomposition theoremZ® is the sum of an infinite 
number of sinusoidal components. Thus the random sea surface elevation may 


be expressed by 


ne H. 
; koa 
= —e a = . + ° (-] 
n ay = sa poit (key - a, t Yao) (3.14) 
AQ+ 0 
where 

ic = an index denoting a component of the frequency spectrum 
Ars = the wave height of component ic 

a the wave number of component ic 


1) 


a random phase angle 
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Let the random phase angles be independent with probability density func- 


tions given by 


LAN 
= 

1A 
| 


ty (vy) = sa 0 
= 0 otherwise (3.15) 


Further, let the wave height of sinusoidal component ic be related deter- 
ministically to the contribution to total energy density at frequency Qe O° 
Pierson!® shows that by the central limit theorem, the sea surface eleva- 
tion which results from these assumptions converges in quadratic mean to a 
Gaussian random process. 

Let bE be the contribution of sinusoidal component ic to the 


total average energy per unit surface area. Then, 


* ne k 
B=) AEs (3.16) 
1Cc= 


The average energy per unit surface area of a sinusoidal wave with height 


H is given by2® 


E = 42 (357) 


* _ Y 2 
aE. = og Hi, (3513) 


Equation 3.13 may also be written 


A nc 
ee lim ~~) Sin. 42 (3.19) 
Nes” jcz] re 


AQ=0 
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From Eqs. 3.16 and 3.19, an alternate expression for the contribution of 


sinusoidal component ic to total average energy per unit surface area 1s 


* 
AES. = Sane AQ (3.20) 


Equations 3.18 and 3.20 yield the following expression for the 
wave height of sinusoidal component ic in terms of the sea surface eleva- 


tion spectral density at Qe Oe 


He. = 2972 Sane AQ (3.21) 


When Eq. 3.21 is substituted into equation 3.14, there results 


nc 
n = J v2 Sn AQ sin (key an van) (3922) 


In the limit as Ag approaches zero and nc approaches infinity, Eq. 3.22 


becomes 
n = [ v 2s d2 sin ( ky - at + ¥) (aec3) 
0 


The next step in the development of the relationship between 


Sint) and Sif (2) is the application of linear wave theory to determine 
nk n& 


an expression for random velocity potential o(x,y,t), which is defined in 


such a way that 


U = - — = the random water particle velocity in the y 


(horizontal) direction 


8x = The random water particle velocity in the x 


(vertical) direction 
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The assumption of incompressible, irrotational flow and continuity re- 


quirements lead to Laplace's equation, 


Vo @ = 0 toca) 


The velocity potential is determined by solving Eq. 3.24 and 


applying the following boundary conditions. At x = 0, 


Z a ng (3.25a) 
nix = h + y 
_ 12 
n= 9 3t (3.25b) 
and 
= (3.25c) 
aX at ‘ 


Equation 3.25a states that the water velocity component perpendicular to 
the ocean bottom is zero at the ocean bottom. Fquation 3,25b, which 1s ob- 
tained from Bernoulli's equation, states that the pressure is equal to 

zero on the free surface. Equation 3.25c is a mathematical statement of 
the condition that particles on the free surface stay there (the vertical 
component of water particle velocity at the free surface equals the time 
rate of change of the surface elevation). In linear wave theory, it 15 
assumed that the wave amplitudes are so small] that satisfying Eqs. 3.25b 
and 3.25c at the mean water level, x = h, is nearly the same as satisfying 


them at xX = h + yn. Thus the last two boundary conditions become 


eer (3. 26a) 


iy = h 


n = 


Q 


1 
g 
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and 
= 2 = 2 (3.26b) 
OX, ulheehe at 
A solution which satisfies Eqs. 3.24 and 3.25a is 
ne 
¢ = bo C.. cosh (k. x) - COS (ks oy - 25 .t + Yeo) (3227) 
Ci is a coefficient which is determined by substituting Eqs. 3.22 and 3.27 


iieorkq, 3.26a. 
C. = 22. ¥ 2S AQ Os ALY (3.28) 


Substitution of Eq. 3.28 into Eq. 3.27 produces the linear wave theory 


solution for the random velocity potential. 


‘ cosh (k; x) 
6 = ee as AQ 
ic=] “ic UM cosh (Kh) 
* COs (key - 25 t + wae (3,29) 


In the limit, as Ag approaches zero and nc approaches infinity, 


[>] 


7 9 V35—Tayda cosh_(kx 
: [ va ot A1en cosh (kh 


6) 


cos (ky - at + ¥) (3.30) 


Substitution of Eqs. 3.22 and 3.29 into Eq. 3.26b leads to the following 


useful relationship. 


—— 
QM. =F YG kK. tanh (k. oh) (aes) 
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The horizontal component of random water particle velocity is 


obtained by differentiating -the velocity potential with respect to y. 


y = 2 28 
eee 


2 ¥ 25), oe sos 7 - sin (ky - at +¥) (3.32) 


oW—8 


The horizontal component of water particle acceleration is the derivative 
of the particle velocity with respect to time. 
ws 


aU : aU 
De ot 7 Moy t MK (3.38) 


The last two terms of Eq. 3.33 are convective acceleration terms which 7 


may be ignored as being small compared to the temporal acceleration, a. 


ce 
i 
aa 
cti/C e 


e rae -~Siccash (kx) 
0 


* cos (ky - at + ¥) (aaea) 


3.3.2 Correlation Functions of Random Water Particle Velocity 


Because the random water particle velocity is a linear function 
of the zero mean, Gaussian, stationary random sea surface elevation, the 
water particle velocity is also a zero mean, Gaussian, stationary random 
process. Therefore, the correlation function of the water velocity at 
nodes nk and n& is independent of time, t, and is a function only of the 


time lag, t. By applying Eq. 3.32, the correlation function may be written 
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FO 
e 
Ce 
4 
a alll 
i 


5 [U,, (t) U p(t + t)] 


cosh (kx) 


2 gk 3s 
; | Q on Gosh (kh) 
- 2=0 


io) 
i 


sin (kY - ant + ¥) | a v AS dn: 


7 -=0 
cosh (k Xp) 


-_ —————— sin (kere -n (t+ +9) | (335) 
cosh (k h) 


The right hand side of Eq. 3.35 contains four random variables - 
va and Ye? the random deflections of nodes nk and n& and ¥ and v, the 
random phase angles. The random node deflections may be eliminated from 
Eq. 3.35 by assuming that the node deflections are small compared to the 
wave length of ocean waves, This assumption is reasonable if the deflections 
are on the order of a few feet, because typical ocean wave lengths are on 
the order of a few hundred feet. Another way of stating the assumption is 
that the deflections are so small compared to the wave lengths that the 
water velocity and acceleration at the undeflected riser node points is 
essentially the same as the water velocity and acceleration at the deflect- 
ed node points. 

By definition of the expectation operator, the correlation func- 


pion of Eq. 3.35 is 
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7 7 - cosh (kx, _) 
es — gk rxe——a9 nk 
RYU ¥} { ij cn >nn ° Cosh (kh) 
N=0 


nk“ n& joo eee 


r cosh (k'x_,) 
a nn cosh (k h) 


- sin (-a(t + 1) + | - fyyt (vab' dy dy (3.36) 


in which Foy! the joint probability density function of the statistically 


independent random phase angles ¥ and y', is 


H 
_s 
-_ 

2 
— 


fuy! (Web ) ; fy: (v) 


= —y7, 0s 0 s 2m, 0O 5 Ww g 2m 


otherwise (337) 


n 
© 
vy 


Equation 3.36 may be rewritten in a more convenient form by interchanging 
the order of integration with respect to phase angle and circular frequency 
and by changing the limits of integration with respect to phase angle so 


as to exclude the region where Foy is zero. When this is done the result 


1S 
r cosh (kx.,.) cosh (k x.,) 
ee k ne 
ee = || iy 2 ark! Vs! da da’ Ye 
(cael a ae nn cosh (Kh) cosh (kh) 
Qu QT 
sin [-at + w] sin [-2 (t+) +0] 
v=0 y =0 
—s dy dy’ (3.38) 
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Consider the double integral with respect to phase angle which 
emgs Eq. 3.38. If a # 2°, the double integral reduces to two single in- 


tegrals each of whose values is zero. Hence, Eq. 3.38 has a nonzero value 


only when Q = 2’ and may be simplified to 
- cosh (kx.,) cosh (kx_,) 
Ry og (cr) = | 2( HK)? 5 (9) ne 
nk ne 4=0 cosh” (kh) 
on Qn 
ay i f sin [-at + vy] sin [-a(t + +) + vl 
tm G0 Y=0 
dw dy da (3.39) 


The double integral with respect to random phase angle, which is easily 
evaluated by writing the product of the two sines as the sum of two co- 


sines, equals oné cos (nt). Therefore, Eq. 3.39 reduces to 


cosh (kx 1) cosh (kx, ») 


By (t) = ‘A (9% s (a) Se 
Unkene 4 a i Me cosh” (kh) 


cos (Qt) da (3.40) 


By virtue of Eq. 3.7, the variance (if nk equals n@) or covari- 
ance (if nk # n£) of the water particle velocity may be obtained by setting 


tT equal to zero in Eq. 3.40. 


Z 
o; e = Rie e (0) 
Vakene Une 
- cosh (kx...) cosh (kx, ) 
| c( ak é Ss (a) ee le J dea (3.41) 
5=0 o cosh (kh) 






By 


3.3.3 Spectral Density of Water Particle Velocity 


The expression for S* *, (2) may be obtained by simply comparing 
Unk une 


Eqs. 3.10 and 3.41. It is apparent that S,, (2) must be the expression 


nk ene 
in brackets in Eq. 3.41. A more rigorous derivation follows. 
By definition, the spectral density function is the Fourier trans- 


form of the correlation function. 


E ] =1QT 
So hae + f mo ee (3.42) 
Unk n& én ; Unk n2 


Substitution of Eq. 3.40 into Eq. 3.42 yields 


[o a) c ioe) 
ae 1 2 i k 2 ! 
S = —_ S 
ies Sez J g _ Soe ye, 
; T==a 2 =0 
cosh eon) cosh (kX) 9) 
5 os(2'r) do 
Gosia tk mM) 
a Ten (3.43) 


in which the prime notation indicates the wave number and frequency terms 
introduced by Eq. 3.40. Interchanging the order of integration in Eq. 3.43 


leads to 


a 2 ,_ cosh (k'x.,) cosh (k x.) 
S50) = Hf (He) 5) (a!) rk nee 


Ae a cosh“(k'h) 


‘| oe Gk ee | do’ (3.44) 
T== 
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Consider the expression in braces in Eq. 3.44. The solution of 


this integral, as given by Davenport and Root?? is 
i cos (a's) e 18? dy = 9 [6 (n= 0') +6 (a -0')] (3,45) 


where 6() is the Dirac delta function, defined as follows. 


b 
ii (x) dx = 1, a<O0O<b 
a 
=U otherwise 
é(x) = 0, ee 
= Oy, x = 0 


It follows that Eq. 3.45 may be restated 


i cos (a't) e 7 dr = og, Q° = £9 
1 (0s otherwise. (3. 46) 


Because of the nature of Eq. 3.46, the integral with respect to 
2’ in Eq. 3.44 is nonzero only when Q' = +9, Therefore, Eq. 3.44 simpli- 


fies to 


~ 


: cosh (kx 1) cosh (kx) 


, (a) = 7(E)% 8 (jal) a 


. (3.47) 
Unk ene 2 nn cosh” (kh) 


and the corresponding one-sided spectrum for the horizontal component of 


water particle velocity is 
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Seer O)) aeons. S° 


(2) 
Une UU 


nk n& 


eosh (kx) cosh (kx_,) 
cig Sl (a), 


cosh* (kh) 


— Oe a =< 0 (3.48) 


3.3.4 Other Fluid Kinematic Spectra 


The spectral densities, S, = (2), S° * (a2), and S$ (Qe 
ene Yak one Unk ene 
could be derived in the same manner as that used for Si U (2), by writing 
nk n2 


the corresponding correlation functions, substituting Eqs. 3.32 and 3.34, 
and performing a Fourier transformation. A shorter derivation is possible, 


however, and will be used here. 


To obtain eat (2) in terms of Sig (2), it is only necessary 
nk né nk n2 
to recognize that 
= Si 9 (2) = 0 (3.49) 
nk ng 


The substitution of Eq. 3.42 into Eq. 3.49 yields 


(1) 


fos U 
be [- bye si eo wit 4. 
: ] i| -12tT 
- iN» Ro eco e dr = 0 (3.50) 
en = Unk ene 


The partial derivative of the correlation function is given by 
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Ri = =. 5 (edie) Ue es 


id 
~ 


et) 
nk ene 


E [U,, (t) Unp(t + 1)J 


si it ec (1) (35) 
Unk ene 


which, when substituted into Eq. 3.50 results in 
5 if Re U (t) e “int dt 
ns nk ng 


© 1 f ° e -i2t 
= 72 s- R (t) e dt (seu2) 
ory Unk ne 
But the two Fourier transforms in Eq. 3.52 are, by definition, spectral 


density functions. Hence, 


Sot (oO) = sn Ss. (2) (3353) 
Unkne Uke ne 


and for the one-sided spectrum, 


Soe yr = 10S en) (3,54) 
Unkune Varene 


In like manner by using the relationships 


3 

—— Ro 7 ile es Rt. ee) (3.55a) 
at UU ng Unk ene 

3 

— Ree (rt) Rae (cy (3.55b) 
” Vine Unkene 


and 
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ae 
— §° ° (fg) = 0 (3,55¢e) 
oT Vane 


it can be shown that 


Sere (Q) sya eS eee | (0) (3.56a) 
Ui ene Uakene 


and 


2 


S (2) = 2 Sy g(a) (3. 56b) 


Unk ene Unk n& 


3.3.5 Modification of the Input Sea Surface Elevation Spectrum by Linear 
Wave Theory 
It has been shown that the fluid kinematic spectral density func- 
tions may be expressed in terms of the sea surface elevation spectral den- 
sity. In a sense, the ocean acts as a linear filter whose input is the sur- 
face elevation spectrum and whose outputs are the fluid kinematic spectra 
at all locations, x, throughout the water column. In this thesis linear 
wave theory is used to model the effect of the ocean on the input spectrum. 
A convenient way to regard the effect of the ocean on the input 
spectrum is to examine the power transfer function, r(2), which is defined 


as the ratio of the output spectral density to the input spectral density.2° 


S (2) 
- _output 
r(9) Siaput 9 (3.57) 


Consider the water velocity and acceleration spectral densities at node X ak 


as outputs with the ocean acting as linear filter. These spectra may be 


expressed as products of an input spectral density and the appropriate power 





| 


transfer function. 


ig oy =) (Xi a= Seen) (3.58a) 


2S RS (ol) ae ame ed camer te) ac Soe Aes) (3.58b) 
Unk U * nk nn 


An examination of Eqs. 3.48 and 3.56b indicates that the two power transfer 


functions are 


kg cosh (kx) 9 


ry (Xb 3 Oy ee 7 ae (3.59a) 
and 
Ry ox C= 9? Tx 2) 
U nk? U nk? 
kg cosh (kx 
ee re Ke (3.59b) 
cosh (kh) 


Figures 3.3 and 3.4 demonstrate the manner in which the ocean 
modifies the input sea surface elevation spectral density. The input spec- 
trum shown is that of Pierson and Moskowitz for a wind velocity of twenty 
knots and is given by Eq. C.4 of Appendix C. The power transfer functions 
and output kinematic spectra shown are for a water depth of six hundred feet 
and for distances below the surface of one, five, and ten percent of the 
total depth, that is, Shee une and sixty feet below the surface. Two 
distinct effects are seen in the figures. First, except for frequencies 
approaching zero, the output kinematic spectral densities decrease as dis- 


tance below the surface increases. Secondly, except for locations very near 





58 


the surface, the ocean filters out the high frequency portion of the input 
Sea surface elevation spectrum. This latter effect becomes more pronounced 
as distance below the surface increases. It may be concluded that the lower 
frequency components of the input spectrum influence the output kinematic 
spectra to greater distances below the surface than the higher frequency 
components do. 

Figure 3.5 shows profiles of root mean square water velocity for 
various wind velocities, a water depth of six hundred feet, and the Pierson- 
Moskowitz sea surface elevation spectrum. Throughout the water column the 
root mean square water velocity, which is the square root of the integral of 
the water particle velocity spectral density function, increases with an 
increase in wind velocity. For a given wind velocity, the root mean square 
water velocity decreases as distance below the surface increases. This ef- 
fect is caused by the reduction in spectral density and the removal of high 
frequency components of the input spectrum by the ocean as represented by 


linear wave theory. 
3.4 Wave Force Spectra 


In this section, the relationship between the fluid kinematic 
Spectra developed in Section 3.3 and the wave force spectra is derived. 
The mode] used to express the random wave force intensity in terms of the 
fluid kinematic parameters is the nondeterministic version of Morrison's 


formula,?> 


Psi) = ste (0) sales U (3.60) 


where 


og 


Ga i ace 5 D (3.614) 
and 
nD° 
C. = Cy e 7 (3.61b) 


In Eq. 3.60, the drag and inertia coefficients, mass density, and diameter 
are deterministic quantities, and the wave force intensity, P(x,t), and 
fluid velocity and acceleration, U and U, are random. 

The random wave force is composed of a drag component which jis a 
nonlinear function of the random water particle velocity and inertia com- 
ponent which is a linear function of the random water particle acceleration. 
The presence of the nonlinear velocity term complicates the solution be- 


cause it precludes the use of a superposition technique. 
3.4.1 The Linearized Wave Force Model 


The troublesome nonlinearity in the velocity term may be removed 
by the method of equivalent linearization, which is described in detail by 
Lin?’ and used by Borgman'’»!® and others!229 to linearize wave drag forces. 
The essence of the method is to replace the nonlinear term, U|U|, with a 
statistically linearized term, C U, where C is a constant. The resulting 


random error, E is 


eve JUN oe FOU (3.62) 


The coefficient C is determined such that the mean square value of the 


error is minimized. 


A ree 
aoe ee) dO (3.63) 
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Because the operations of differentiation and determining the expected 


value are both linear, Eq. 3.63 may be written 
pie =] = 0 (3.64) 
aC i 
The substitution of Eq. 3.62 into Eq. 3.64 yields 
°? e 2? 7 
ae Ee leet ee eat)  =950 (3.65) 


Because U is a zero mean random variable, E a = oF , the variance of the 


velocity. Therefore Eq. 3.65 may be written 


2 J 
Gs 5 (3.66) 


The numerator of Eq. 3.66 is evaluated by noting that for U < 0, 
\U| = - U, and for U > 0, |U| = U. Hence, 


f U7 jul fy (a) du 


MW 


E CUe|U | J 


+ i of e(uledu (3.67) 


J 


where fj, (u) is the normal probability density function given by 
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Substitution of Eq. 3.68 into Eq. 3.67 and evaluation of the integrals leads 
to 


E(U2jU) = VJ of? (3.69) 


When Eq. 3.69 1s substituted into Eq. 3.66, there results 


8 
g (3.70) 


tl 


The value of C determined by Eq. 3.70 minimizes the mean square random error, 


E, because 
3° 2 oC; o£ 
= 2€ [U2] 
2 ne 
=e Of, (3271) 


and the variance, of » must be nonnegative. The linearized version of the 


wave force equation may now be written 
P(x,t) = AVE OF) U(x,t) + Cc. Ui) (272) 


Because the random force intensity given by Eq. 3.72 is a linear function of 
the zero mean, Gaussian, stationary velocity and acceleration, the force in- 


tensity 1s also zero mean, Gaussian, and stationary. 
3.4.2 Correlation Function of the Random Wave Force Intensity 


The correlation function of the random wave force intensity at node 
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nk, Pact)» and the random wave force intensity at node nf, Papit)s is by 


definition 
Pp Ge esl Pal Ge): Foy (t + 1)] C3573) 
where 


P K (ti =. 4P (x) ot) 


Substicucion of Eq. 3.72 into Eq. 3.73 leads to 


De) 
3 [co 


Z 
ee Ge Re es (t) (3.7 4) 
Unk Une 


3.4.3 Wave Force Spectral Densities 


The wave force spectral density is obtained by taking the Fourier 


transform of the corresponding correlation function. 


tel 


] -i2t 
S (Q) = = { R (2) e dt (3375) 
Pak ng ae es a 


Equation 3.75 may be expanded by substituting Eq. 3.74. 
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] 8 
(2) = = i 1G = ee con Co, Rae sic (+c) 
PoP on uy Ue a Unk Une Une 


Pe Se Pe ea (3.76) 
Unk ene 


Interchanging the order of summation and integration in Eq. 3.76 leads to 


2 8 
5 = Cpa a ery Scene 2) 
ae ne ai Unk n& Unk nl 
+ CC Oo, See. a) 
ue Unk Vik Yne 


eS) 
a Unk ne 


Equation 3.77 may be further simplified by employing Eqs. 3.54 and 3.56. 


2°68 ; 8 
S (oy = C9 = =m econ. “Teeee \/2 (co; 67) 
Pak ng uot Unk Une Bneey Uae ~ “Une 
2 2 
+O Se a) (3.78) 
: Vane 


Equations 3.78 and 3.48 constitute a mathematical link between the input 
to the problem - the surface elevation spectrum - and the spectral density 


function of the wave force intensities. 
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3.4.4 Introduction of the Lumped-Smeared Force Model 


The accuracy of response calculations may be improved by incor- 
porating the lumped-smeared force representation developed in Chapter 2 
into the nondeterministic force model. Consider Eq. 3.78, which with the 


use of Eq. 3.48 may be written 


+ & cyl Mis (a) 


cosh (kx, ) cosh (kx, >) ean 
cosh (kh) cosh (kh) 

The last two terms in Eq. 3.79, the fractions involving the hyperbolic 
cosines, represent the spatial variation of the random force intensity 
along the riser axis. The term cosh tee is also a spatial modifier 
for the inertia force, as indicated by Eq. 2.26. In the lumped-smeared 
force system, the ratio of hyperbolic cosines is replaced by the spatial 
modifier, Any (2) » which is given by Eq. 2.35. The lumped-smeared force 
representation is introduced into the nondeterministic model by replacing 
the hyperbolic cosine ratios in Eq. 3.79 with the lumped-smeared spatial 


modi fier, An (2) » of Eq. 2.35 resulting in 
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2 8 
S OQ = Cee. Ew ce on 
Pak ng T Unk Une 

: 8 Ze 

aoe) (Go ye (o = 07 )) + C.- 20am 

ua 1 UK Une a 

Pot i icemiaeceein) s.. (a) (3.80) 

Q Gnk Gne nn ; 


3.5 Response Spectra 


The final step in the calculation of the nondeterministic static 
response of the riser is that relating the random wave forces to the random 
response. Response parameters of interest are the node deflections and mo- 
ments and the bottom rotation of the riser. The approach adopted here is to 
first calculate the spectral densities of the node deflections and express 


the remaining response spectra in terms of these deflection spectra. 
3.5.1 Random Node Deflections 


Let {Y(t)} be an n by 1 vector of random node deflections caused 
by the random wave forces, {P(t)}. The random deflections are related to 


the random force intensities by 
CK] {Y¥(t)} = {P(t)} (3.81) 


in which the stiffness matrix [K] is defined by Eqs. 2.12. To obtain {Y(t)}, 
Eq. 3.81 is premultiplied by [kT!, 


y(t)? = (KT! cece) (3. 82) 


A form of Eq. 3.82 which is useful for calculating the correlation functions 
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and spectral densities of the deflections is 


n 
Tg (eee) Keer) (3.83) 
k 


Because the node deflections are linear functions of the random force in- 
tensities, they are zero mean, Gaussian, and stationary. 
By definition, the correlation function of the response at nodes 


ni and nj is given by 


R Se et at Catt ea (3.84) 


eal nJ 


ni nj 
When Eq. 3.83 is substituted into Eq. 3.84 and the order of summation and 


expectation is interchanged, there results 


It =| = 
y L K K R (tr) (3,85) 


n 
Gey Sa : 
vf k ni,nk ‘nj ,nd Parone 


ni nj nk=1 n 


The response spectral density is obtained by a Fourier transfor- 


mation of the correlation function. 


sf ] -1QT 
g (ay eee fr ee ‘e (3,86) 
eae a ae 


Substituting Eq. 3.85 into Eq. 3.86 and interchanging the order of integra- 


tion and summation leads to 


n n 
: 2] =| 
oa ) ) K.. Ko S (2) (e827) 
‘nt 'ng AE RES ni,nk ‘nj,né BAL 
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Equation 3.87, together with Eq. 3.80 may be used to calculate the spectral 
density of the random static deflection resulting from the random sea des- 


cribed by the sea surface elevation spectrum, oe (oie 
3.5.2 Random Bottom Rotation 


Consider the random bottom rotation, O,(t) defined as the slope 


of the deflection curve, Y(x,t), at x = 0. 


a,(t) = & Y (0,t) (3.88) 


By using Eq. 2.3, a finite difference version of Eq. 3.88 may be written, 


a(t) = Sa Lv, (t) + Y(t) J (3.89) 





where vy is the deflection of an imaginary node located a distance z : j 


below the riser bottom. Because the bending moment at the bottom of the 


riser is zero, Y_, (t) = SF (t), as shown by Eq. 2.9b, and 
_ nt] 
oS a Y, (t) (3.90) 


The autocorrelation function of the bottom rotation is 


i (7) = El ait) ome (ieateam (3:91) 


which, after substitution of Eq. 3.90, becomes 


(rt) (3.92) 





R ee: 
(y= (RELY Ry 


A Fourier transformation of Eq. 3.92 yields 
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(2) (3.93) 





3.5.3 Random Bending Moments 


Like the bottom rotation spectral density, the bending moment spec- 
tral density may be determined from the deflection spectral densities. Let 
M(x,t) be the wave-induced random bending moment at coordinate x and M4 (t) 
be the random bending moment at node ni. From small deflection beam theory, 


the bending moment of the riser is related to the deflected shape, Y(x,t), by 


M(x,t) = - El ae (3.94) 


Equation 3.94 is transformed into finite difference form by introducing 


Eq. 2.2 for the second derivative. 


- n+1 2 
M at) = at ( L : ) EI [ Vey = ay * + i 


ni ni nit+l J (3.95) 


This transformation between random node deflections and random node moments 


1s conveniently expressed in matrix form, 
ee = Cry} (3.96) 


where [ J ] 1s an n by n, three term banded transformation matrix whose 


nonzero elements are 











Ing ni = 2EL( Pe ean en (3.97a) 
_ ese les? ae 

are = - EI ( 1 eee eet (3.97b) 

Innie) 7 7 EL (SEED, oni = aan (3.97c) 
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A convenient form of Eq. 3.96, which takes account of the banded 


nature of [ J ], is 


nit] 


Mo eye ee a ent en (ct) (3.98) 


The correlation function of the random moment at node ni and the 


random moment at node nj is 


Maing (tr) = E [Ms (t) Mag (tet a) (3.99) 


Substitution of Eq. 3.98 into Eq. 3.99 leads to 


i ee (1) (3,100) 
R (tT) = ae ese R T 3. 100 
Mi nj nk=ni-1  n@=nj-1 (M7 enk ndsne Viet ne 


A Fourier transformation of both sides of Eq. 3.100 results in 


nit] ng 
) Pie tlhe 9: MRIS (2) (3.101) 
nj nk=ni-1  ng=nj-1 M19MK ng ane “YAY np 


3.6 Sample Problem 


In this chapter, a model for calculating the response of a marine 
riser to random wave forces has been developed without consideration of 
dynamic effects. In order to illustrate the influence of nondynamic effects 
on the response, a sample problem has been solved with the use of this 
model. 

The riser installation chosen for the sample problem has the 


following characteristics, 
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G =so2 10 

D = 20 inches 

h=L = 609 feet 

E = 29 x 10° pounds per square inch 
1 = .0621 feet* 

W = 248 pounds per foot 


The ratio of top tension to riser submerged weight is 1.2, and drag and 
inertia coefficients are 1.0 and 1.5 respectively. A finite difference 
model with 31 equally spaced nodes is used to represent the riser, and the 
random sea surface elevation is taken as that described by the Pierson- 
Moskowitz spectrum. 

If the entire mathematical model is considered to be a linear 
system whose input is the sea surface elevation spectral density function 
and whose output is a response spectrum such as the bottom rotation spectral 


density, then the problem may be expressed as 


S (a), = Ee Ag)e=eS oc a(n) (37102) 
where ly (2) is the bottom rotation power transfer function of the system. 
0 
From Eqs. 3.80, 3.87, and 3.93, it follows that the bottom rotation power 


transfer function is given by 


7) 


5 i Q 
n n 
-| -| 2.8 

: \ K K We Or ae ae 

tet noe) Lenk alieme uot Une 
cai) cic s (oj - Of ) 

nk nz 

+ 2 92 ] q (2) a., (a) (3.103) 

a Qnk She 


Figures 3.6a, 3.6b, 3.6c and 3.6d show the input sea surface ele- 
vation spectral density, the bottom rotation power transfer function, and 
the static bottom rotation spectral density for wind velocities of 10, 20, 
30, and 40 knots. The variation of variance and standard deviation of the 
bottom rotation with wind velocity is shown in Fig. 3.7a and 3./7b. 

In all cases, the bottom transfer function decreases markedly with 
increasing frequency. It may be concluded then that the effect of linear 
wave theory, the Morrison force formula, and the structural response to 
static wave loads, is to filter out the high frequency components of the 
sea surface elevation spectrum. This is illustrated in Fig. 3.6a through 
3.6d where, in each case, the output bottom rotation spectrum has a narrower 
frequency band (is more sharply peaked) than the input sea surface elevation 
spectrum. 

For low wind velocities, the static response of the riser jis al- 
most entirely a result of the inertia force. At a wind velocity of 30 knots, 
that part of the response caused by drag force is perceptible, and at a wind 
velocity of 40 knots, the effect of the drag force is appreciable. The ra- 
pid increase in the importance of the drag force at higher wind velocities 


is a Consequence of its nonlinear nature. 
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Chapter 4 
THE DYNAMIC RESPONSE OF A MARINE RISER TO RANDOM WAVE FORCES 
4.1 General 


A major objective of this thesis is the development of a mathe- 
matical model for predicting the random dynamic response of a marine drill- 
ing riser to the random forces produced by ocean waves. Parts of this model 
are derived in the preceding two chapters. A finite difference structural 
model and a lumped-smeared force model are developed in Chapter 2 and are 
combined with linear wave theory and Morrison's force equation in Chapter 3 
to formulate a model for determining the static response of the riser to 
random sea waves. Dynamic considerations are purposely avoided in Chapter 3 
in order to simplify the overall! development and focus attention on how the 
input sea surface elevation spectrum is modified or filtered by linear wave 
theory and Morrison's force formula. In this chapter, response spectral 
density functions are derived with dynamic effects considered. 

Some of the relationships derived in Chapters 2 and 3 for the 
Static riser problem remain valid for the dynamic problem. The power trans- 
fer functions relating fluid kinematic spectra to sea surface elevation spec- 
tra are equally applicable whether the problem is static or dynamic. As 
Shown in Section 2.4.2, the natural frequencies and mode shapes of the riser, 
which are needed in the dynamic analysis, may be obtained by using the finite 
difference model. The lumped-smeared system of representing wave force in- 
tensity 1s tndependent of time effects and thus is also applicable to the 


dynamic problem. 
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A unique feature of the dynamic riser problem is that the random 
wave force intensity is dependent upon the random riser response. The 
hydrodynamic force intensity 1s a function not of the absolute water par- 
ticle velocity and acceleration, but of the relative velocity and accelera- 
tion of the water particles with respect to the structure. In addition to 
modifying the wave force, wave structure interaction also produces a hy- 
draulic damping effect on the structural response. 

The method used here for determining the dynamic response of the 
riser to water waves is adapted from recently published works of Foster!? 
and Malhotra and Penzien,29;2! who developed methods for calculating the 
random response of offshore towers to ocean waves. The method is modified 


as necessary to meet the requirements of the riser problem. 
4.2 Equation of Motion for Dynamic Riser Problem 


The governing differential equation of the marine riser is Eq. 1.3, 


which is repeated here for convenience 


4 2 
eI SY - 1(x) SY - w = p(x) (1.3) 
dx dn dx 


When dynamic effects are significant, the riser deflection, y, 1S a function 
not only of position, x, but also of time, t, and the derivatives in Eq. 1.3 
are properly written as partial derivatives. 
3 3¢ 3 
EI <p ASE), hie) ys mee = y(x,t) 5 p(x,t) (4.1) 
5 4 2 OX 
X OX 
The lateral force on the riser, p(x,t), may consist of several 


components-~a random hydrodynamic force, a force whose effect is equivalent 


to that of a random top offset, a d'Alembert force caused by opposition of 
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the riser mass to acceleration, and a damping force which is also opposed 


to the motion of the riser. 


pct) = BCC rm y. Py(xst) (4.2) 
where 
P(x,t) = the random hydrodynamic force 
P(x,t) = the force intensity equivalent to a random top 
offset 
m = riser system mass per unit length 


Py(x,t) = the damping force 
4.2.1 Equation of Motion in Matrix Form 


When the method of finite differences is applied to Eq. 4.1, the 
riser is replaced by a system having n equally spaced nodes and n degrees 
of freedom. The terms on the left side of Eq. 4.1 are gathered into the 


stiffness matrix defined by Eqs. 2.12, and Eq. 4.1 assumes the form 


ia ee (V(t) = (pCa (4.3) 


nxn 


In Eq. 4.3, the force vector {p(t)} may be replaced by its constituent 


forces as given in Eq. 4.2. 


ip(t)} = 4P(t)} = aP(eiy = | mel oli: (4.4) 
where 

{P(t)} = a vector of random hydrodynamic force intensities 

{P(t)} = a vector of force intensities which are equivalent 


to the top offset 
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[--mJ = an nxn diagonal mass matrix 

el = an nxn matrix of damping coefficients 

{Y} = a vector of random riser accelerations at the nodes 
{Y} = a vector of random riser velocities at the nodes 
[c]t¥} = a vector of damping force intensities 


Consideration of the effect of random top offset is deferred until Chapter 5. 
For the present, {P{t)} is taken as zero. 

The applied random hydrodynamic force intensity is a function of 
the fluid velocity and fluid acceleration, as indicated by Eq. 2.23. The 
fluid velocity may consist of a steady current and an unsteady component 
caused by surface waves. The effect of a steady current is studied in Chap- 
ter 5. For the present discussion, the steady current is taken as zero, and 
the fluid velocity 1s solely that caused by waves. If the structure responds 
to the applied force in such a way that the structural velocities and acceler- 
ations are significant with respect to the fluid velocities and accelerations, 
then the wave forces of Eq. 2.23 should be written as functions of the rela- 
tive fluid velocity and acceleration with respect to the structure. 

Let V be the relative velocity of the water with respect to the 
structure and V be the relative acceleration of the water with respect to 


the structure. 
a ie (4.5a) 
Sl) 3 (4.5b) 


The variables V and V are both linear functions of zero-mean Gaussian random 


variables and are therefore also zero-mean Gaussian random variables. 
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When wave structure interaction is considered, Morrison's formula, 


Eq. 2.23, assumes the form 


2 
nee Dory , 
P(x,t) = Cp eDV{V| + Co a= V (4.6) 


or in matrix form 


(P(t)} = Ca{V|V|} + C.1V} (4.7) 
where 

Uda ~ (Unk 7 Uae er 7 ara 
and 

‘nk : Unk 7 Vnk 


Equations 4.3, 4.4, and 4.7 combine to yield the following 
[-mJi¥} + [c]t¥} + [K](Y) = C.CV|V|} + C03 (4.8) 


Equation 4.8 represents, in matrix form, the n simultaneous differential 
equations of motion for an n degree-of-freedom finite difference model of 

a marine riser subjected to random wave forces. The drag force term is a 
nonlinear function of the relative velocity. This nonlinearity will be re- 
moved by the method of equivalent linearization which was introduced in 


Section 3.4.1. 
4.2.2 Linearization of the Equation of Motion 


The technique of linearization used here is essentially the same 
as that used by Malhotra and Penzien.2° First the structural displacement, 
velocity, and acceleration terms of Eq. 4.8 are removed by the trans forma- 


tions 


Ue 


{Y} = {U} - {V} (4.9a) 
epee Se an (4.9b) 
{Y} = {U} - {V} (4.9c) 


In Eq. 4.9a, {U} and {V} are vectors of the absolute and relative fluid dis- 
placements at the nodes. Equations 4.9 follow from Eqs. 4.5. When Eqs. 4.9 
are substituted into Eq. 4.8 and the resulting terms are rearranged, the n 
simultaneous differential equations are expressed in a relative coordinate 


system with respect to the riser structure. 
Pm ](V} + Cove + CK](V} + CVV] = Pmt} 
+ [c](U} + [K]{U} (4.10) 


where 


fees) =) Priel ee te 


and [~I~] is the unit diagonal matrix. 
In Eq. 4.10, there is one nonlinear term involving the relative 
fluid velocity with respect to the structure. This nonlinearity will be 


removed by the method of equivalent linearization. Let 
[c]iv} + ccv[V]} = [e]tv} + ce} (4.11) 


where [c] is a matrix of modified damping coefficients and {&} is a vector 
of the errors which result when [c iV} is used to replace the two terms on 


the left side of Eq. 4.11. 


(E} = (Lel - [c])(V} + C-eviv|} (4.12) 
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Elements of the error vector are functions of the random variable V, and 


hence are random variables themselves. It 1s convenient to let 


anna = oe “a for alee ny (4.13) 


that is, the off-diagonal terms of [c] are assumed to be identical to the 
off-diagonal terms of [c]. This is a reasonable assumption because it 
simply implies that the hydraulic damping force at a given node is inde- 
pendent of the relative velocity at other nodes. This assumption permits 


an uncoupling of Eq. 4.12, which may now be written 


; . 7 ‘ = Sey (G0 5 
Cniond ~ Snjanj) “ng * CaM ng!Y 


= ( | (4.14) 


Ej nj 
To linearize Eq. 4.10, elements of the modified damping matrix 
[c] are selected such that the mean square values of the elements of the 


random error vector, {E}, are minimized. The mean square error is the second 


moment of the error or the expected value of the square of the error. 


2 _ 0 4 ° 2 
ELE, «J = E[{(cy; 3 - Ciena nd + Gna onae 1] (4.15) 
The mean square error iS minimum with respect to on a when 
= (ELEC s}} oa (4.16) 
~nj nj 


In Eq. 4.16, the partial derivative and the expectation, both being linear 


operations, may be interchanged. 


9E_. 
EE, spot! = 0 (4.17) 
J <nj nj 


(2 


From Eq. 4.14, 


aE. : 
a= See (4.18) 
Thiet) J 


Substitution of Eqs. 4.14 and 4.18 into Eq. 4.17 results in 


Belle J vat Cav [11-513 = 0 (4.19) 


= Cn. ; 
NJ »NJ “NIN” NAY U nj nj 


Equation 4.19 may be rearranged by combining terms and taking constants 


outside the expectation operator.. 


yd a7 e , 
Brana’ scndenae © Cong aa CyELV AS Yng Hd = 0 (4.20) 
Solution of Eq. 4.20 for Sar at results in 
©? ° 
Slee | Woe |i 
= Wee li | 
Snjonj (4.21) 


opeprmmenerttes 960 
nen u a2 
ELV, 54 


Because the relative fluid velocity at node nj, Vinge? 1s a zero- 


mean Gaussian random variable, 


2 _ 2 
ELV, 5] - ihe (4.22) 


From Eq. 3.69, it follows that 


a ee ie: 
ELV Wagl] = NEE: . (4.23) 


nJ 


Substitution of Eqs. 4.22 and 4.23 into Eq. 4.21 leaves 


=o | 
“nj sng Cnjsng * UVa ae (4.24) 
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Equations 4.13 and 4 24 are conveniently expressed in matrix form. 


el = Cel +o,\/2 Poy (4.25) 


where [oy J 1S an nxn diagonal matrix of standard deviations of the rela- 
tive fluid structure velocity at the riser nodes. 
The modified damping matrix, [c¢], given by Eq. 4.25, minimizes 


the mean square error because 


nj 


2 see 
ao E[eea = £ eae = E[V*.] = of | (4.26) 
Cnn ty 


OC ea. 
“nj sn 


and the variance 1s always non-negative. With the mean square error mini- 


mized, an approximation to Eq. 4.11 may be made by dropping the error vector. 
[c]iV} + C- (VIVJ} > Cet) (4.27) 


Equation 4.27 1s suitable for substitution into the matrix form of the 


equations of motion, Eq. 4.10. 
Cm1iV} + [clt¥i + [kK]tV} = PmJtu} + [c](0) 
+ [K]iU} (4.28) 


Equation 4.28 may be transformed from the relative coordinate sys- 


tem back into the absolute system by utilizing Eqs. 4.5 and the relationship 
WN SU eT (4.29) 


which follows from Eq. 4.9a- This transformation leads to 


8] 


Pom, JeV} + [e]iv} + CK}CY} = (Dm J - bmJ)tU} 


+ ([c] - [c]) (U3 (4.30) 
But 
Dm] - tn = ¢.t1 
and 
ied) — 6 VO tea 
Therefore, 


ee '. | q oe | ; Bap see ° 
[~m, ~]{¥} + felis + (KY) a C {Us + mate oy 1{U} (4.31) 


Equation 4.31 is the linearized version of the nonlinear equation of motion, 
Eq. 4.8, relating random riser node displacements, {Y(t)} to random water 
particle velocities and accelerations, {U(t)} and {U(t)}. 

In Eq. 4.31, the drag force and the modified damping matrix are 
functions of the matrix of standard deviations of the relative fluid velo- 
city, which is itself a function of the response. Therefore, an iterative 


procedure is required for the solution of Eq. 4.31. 
4.2.3 Technique of Solution of Linearized Equations of Motion. 


The equation of motion in matrix form, Eq. 4.31, actually repre- 
sents n simultaneous differential equations. To obtain a solution relating 
{Y} to the random water particle kinematics, a choice of techniques is 
available. Solving a similar problem for an offshore platform model with 


five degrees of freedom, Foster!? employs the method of nonclassical, 
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complex-valued mode superposition as outlined by Hurty and Rubinstein.” A 
more convenient and efficient method for the riser, which has many more deg- 
rees of freedom, is classical mode superposition. With a technique similar 
to that of Malhotra and Penzien,“° the n equations of motion, Eq. 4.31, may 
be uncoupled by using the orthogonal properties of the classical normal modes 
Uncoupling the damping terms requires some additional approximations to the 
damping matrix. 

Use of the normal mode superposition technique requires know! edge 
of the structure's natural frequencies of vibration (eigenvalues) and cor- 
responding mode shapes (eigenvectors). In Section 4.3, solution of the 
eigenvalue problem is discussed. The following section (Section 4.4) then 
deals with the application of the eigenvectors and eigenvalues to the pro- 


blem of uncoupling and solving the equations of motion. 
4.3 The Riser Eigenvalue Problem - Natural Frequencies and Mode Shapes 


To determine the natural frequencies and mode shapes of the marine 
riser, 1t is necessary to solve the undamped free vibration problem, which 
is stated mathematically by writing Eq. 4.31 without the damping term and 


with no applied force. 

[m~Jiv} + (K]{Y} = {0} (4. 32) 
it 1S assumed that Eq. 4.32 possesses solutions of the form 

a = Acos (wt + wv) (4.33; 


The double differentiation of Eq. 4.33 with respect to time produces the 


useful relationship 
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: 2 2 
ae =  -w wend (4.34) 


Upon the introduction of Eq. 4.34, Eq. 4.32 becomes 

-Pm] w®(¥} + [K]{Y} = {0} (4.35) 
When Eq. 4.35 is premultiplied by Sok, there results 

Ww 
“kT Cmdr + 0 = 10 (4.36) 
W) 

Let [D] = [K]”'C'm.] be the dynamical matrix. Then, Eq. 4.36 simplifies 
to the characteristic equation 

[DIYs = 4 «Y} (4.37) 

U) 

Solutions to Eq. 4.37 are the natural frequencies, Wee and corresponding 
mode shapes, cD ae of the marine riser. 

Eigenvalues and eigenvectors of the marine riser were calculated 
with the University of I]linois subroutine EIGENP, described in Chapter 2. 
The favorable agreement of the eigenvalues for the first three modes with 
other published solutions is also discussed in Chapter 2. 


The solution for the neh mode satisfies Eq. 4.37 


coyie'")} = Legh); 
Yr 


S 


Uy tm ee) = Light); (4.38) 
r 


= 


Premultiplication of Eq. 4.38 by we [K] leaves 
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we tm, Jeo!" ime: (4.39) 


Equation 4.39 will prove useful in uncoupling the equations of motion for 


forced vibration in Section 4.4. 
4.3.1 Behavior of Riser Frequencies 


It is instructive to compare the fundamental frequencies of marine 
risers to the fundamental frequencies of other structures which belong to 
the same family. Such other structures include the constant tension beam, 
the two limiting cases of which are the simple beam and the unstiffened 
(constant tension) string, and the hanging chain (or string with variable 
tension), which is the limiting case of the riser when the flexural stiff- 
ness approaches zero. Figures 4.1 and 4.2 show the variation of the funda- 
mental frequencies of these structures with the cube root of the length 
ratio, G for top tension ratios of 1.0 and 1.2 respectively. For the 
constant tension beam and the unstiffened string, corresponding average 
tension ratios of 0.5 and 0.7 were used to compute the natural frequencies 

From Figs, 4.1 and 4.2, it 1s possible to make several observa- 
tions regarding the behavior of the fundamental riser frequency. The con- 
stant tension beam frequency constitutes an upper bound on the fundamental 
riser frequency. A lower bound is given by the larger of the hanging chain 
frequency or the simple beam frequency. As the tension ratio increases, 
the difference between upper and lower bounds decreases. As G. decreases, 
the fundamental riser frequency approaches that of the constant tension 
beam; as G. increases, the fundamental riser frequency approaches that of 
the hanging chain (A conclusion similar to this latter observation was 


reached by Frohrib and Plunkett® who found that bending stiffness has a 
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small effect on the natural frequencies of long drill strings.) These 
observations are useful in estimating the fundamental riser frequency if 
that frequency is unknown. 

It is also of interest to study the variation of the first few 
riser frequencies with water depth for a given riser section. Figures 
4.3 and 4.4 show this variation for the first five frequencies of a typical 
riser section and for top tension ratios of 1.0 and 1.2 respectively. In 
general, the frequency for a given mode decreases as water depth increases 
and as top tension ratio decreases. As water depth increases, the differ- 
ence between the natural frequencies of various modes decreases, that jis, 
the natural frequencies become more closely grouped along the frequency 
axis. The significance of this latter observation is that for a given riser 
section and a given sea surface elevation spectrum, the number of modes which 
make an appreciable contribution to the dynamic response tends to increase 


as water depth increases. 
4.3.2 Natural Frequencies of the Riser With Static Curvature 


The preceding discussion regarding riser frequencies applies to 
risers with no static curvature. In practice, real risers almost always 
have some static curvature because of the top offset which results from 
lateral motion of the floating surface support platform. Frohrib and 
Plunkett® have shown that a drill string suspended in such a way that it 
deforms under its own weight has natural frequencies which differ from those 
of a vertical drill string. It will be shown here, however, that for top 
offsets which are consistent with small allowable bottom rotations, the 
natural frequencies of the riser may be taken as those for the case of no 


static curvature. 
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According to Frohrib and Plunkett,® the characteristic equation 
for the drill string with static curvature and with vibration limited to 


the plane of that static curvature is given by 





d"y _d’y dTdy_ 2 d°y. | 
dx dx dx 


where y(x) is the dynamic component of motion and i is the static deflection. 

Because the axial tension, T, is not necessarily vertical, dT/dx is not 

exactly equal to w, the submerged weight of the string per unit length. 
Consider the axial tension in a marine riser. The maximum slope 

in a riser with top offset occurs at the lower end and is limited to a few 

degrees by the nature of the bottom joint. Even if the bottom rotation is 

as much as five degrees, the vertical component of the tension, i is nearly 


equal to the total tension, T. 


US = "1 Cos. (57) ) =" 0.996 1 -2aer (4.41a) 
Furthermore, 
dT 
qi ] oe . 


Thus, if the bottom rotation is limited to a small angle, the axial tension 
may be considered as being essentially vertical in direction and the deri- 
vative of axial tension with respect to the vertical coordinate may be taken 
as the riser submerged weight per unit length. 

The last term in Eq. 4.40 is a function of the static curvature, 
dy /dx®. The magnitude of the static curvature for a given top offset 


depends, to a great extent, upon the length ratio, Ge = wL?/EI. If the 
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riser 1S infinitely stiff, then G. ts zero, and a top offset results ina 
rigid body rotation with constant slope along the riser and no curvature. 
As stiffness decreases and G, increases, a given top offset results in an 
increasingly large static curvature, and the last term of Eq. 4.40 assumes 
greater importance. There must be some range of values of Gy commencing at 
zero and extending to some upper limit, in which the curvature term has a 
negligible effect on the natural frequencies. If the riser length ratios 
are within this range, then the solution to the eigenvalue problem for the 
vertical riser is a very good approximation to the solution for the riser 
with static curvature. 

Another approach may be used to demonstrate the applicability of 
the vertical riser solution to the static curvature case. In their investi- 
gations, Frohrib and Plunkett® found that even for unstiffened strings, if 
the ratio of the horizontal reaction to the weight of the string is small 
enough, then static curvature has little influence on the natural frequencies. 
Figure 4.5 shows the variation of this horizontal reaction ratio with bottom 
rotation for some typical risers with static curvature caused by a top off- 
set. For a given rotation, the horizontal reaction ratio decreases as stiff- 
ness decreases (that is, as Ge increases). The horizontal reactions corres- 
ponding to riser rotations of a few degrees are in the range where Frohrib 
and Plunkett found that static curvature has little effect on the natura! 
frequencies. 

It may be concluded therefore, that if the riser bottom rotation 
is limited to a few degrees, then the eigenvalues and eigenvectors for the 
vertical riser are close approximations to the eigenvalues and eigenvectors 


Of the riser with top offset. The analysis of the riser with random top 
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offset 1s greatly simplified by this situation. It should be remembered, 
however, that for a riser installation in which large bottom rotations are 
permitted and do occur, the use of the vertical riser eigenvalues and eigen- 


vectors for analysis may result in unacceptable errors. 
4.4 Solution of the Equations of Motion 


The solution of the linearized equations of motion, which are 
represented in matrix form by Eq. 4.31, 1s accomplished in two steps. 
First, because the n equations are simultaneous, it is necessary to perform 
a transformation which uncouples them. The n independent equations which 
result from the transformation may then be solved using conventional methods 


of vibration analysis. 
4.4.1 Uncoupling the Equations of Motion 


It 1S convenient to replace the hydrodynamic force terms in 


Eq. 4.31 with a single random force intensity vector, {P(t)}. 
eae CU cal ® [oy J{U(t)} (4.42) 
Let 


[o] > Cte); pg(2, 2 gg ht yy (4.43) 


nxnm 


be a matrix of node deflections for all modes beginning with the fundamental 
mode and including the mode nm, where nm is the number of modes considered 
In the problem. The total number of modes for the finite difference model 
of the riser structure is n, so that nm can take any value from 1 to n. In 


many problems, only the lower modes make a significant contribution to the 
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structural response, and for that reason nm is chosen to be something 
less than n. 
A useful property of [¢] is its orthogonality with respect to 


the effective mass matrix. 
Tire arg 
[ol DmJle] = DW (4.44) 


where to)! is the transpose of matrix [¢]. For the risers studied herein, 
the effective mass per unit length is a constant over the length of the 
riser. Furthermore, the EIGENP routine used to determine the riser eigen- 


vectors produces an eigenvector matrix which is scaled in such a way that 
il - 
(el ted), = tet) (4.45) 


These two conditions permit a simplification of Eq. 4.44. 


Col omaIte] = m Col’ S 106] 


th 


ma [o]'Lo] 


Ma [rte] (4.46) 


Equation 4.46 will be employed later in this section. 
In the normal mode superposition method, the node deflections are 


expressed as a combination of contributions from nm normal modes. 


OR Ge iL {Z(t)} (4.47) 


nxnm nmx1 


where Z(t) is the normal coordinate for the rth 


mode of vibration. The 
variable Z(t) is a random variable; the normal mode matrix, [o], is not 


random, Equation 4.31 may be transformed into the normal coordinate system 
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by means of Eq. 4.47. 
[m,~I[el{Z} + [ello}iZ} + CKI[e]{Z} = 4P) (4.48) 


The useful relationship expressed by Eq 4.39, which is true for 


all modes, r, ranging from | to n, may be written in matrix form as 
[K]fe] = Cmte] Det (4.49) 


where Peed 1s a diagonal matrix of the squares of the natural circular fre- 


Z th 


quencies and Wer 1s the square of the r~ natural frequency. The position 


of we] in Eq. 4.49 is critical. Column r of [o] represents the rth eigen- 
(r)) ( 


vector {¢ For each element of i¢ r)y to be multiplied by we, it is nec- 


essary that Pu] follow [¢] on the right hand side of Eq. 4.49. Substitu- 


tion of Eq. 4.49 into Eq. 4.48 eliminates the stiffness matrix. 
[m,Jtel{Z} + Cellolt2) + Cm eIbu tz} =P) (4.50) 


The key to uncoupling Eq. 4.50 is to premultiply each term of the 
equation by Tae 


Col! EmgJLe1Z} + Co]'Le]loltZ) 
+ fel Dm dtelbotz? = oI! tr) (4.51) 
By means of Eq. 4.46, Eq 4.51 may be simplified to 
m,iZ} + [C iz} +m, wiz} = [oe] 'sP (4.52) 
where 


ic ee le! [cltel (4.53) 


0 


9] 


is the modal damping matrix. All of the terms of Eq. 4.52 are uncoupled 

except the damping term. Unless [CJ 1S a diagonal matrix, Eq. 4.52 remains 

a system of nm simultaneous differential equations and all that has been 

achieved so far 1s to reduce the number of unknowns from n to nm. The 

modal damping matrix, [CJ]: is not, in general, a diagonal matrix, although, 

as it is shown in Section 4.6, the diagonal terms will predominate over 

the off-diagonal elements when damping is treated in a suitable manner. 
Malhotra and Penzien2® present a useful technique by which the 

modal damping matrix may be diagonalized to complete the uncoupling process. 


Let 
[C.1iz} = Sie ieee: We (4.54) 


where Ec .] 1S a diagonal matrix of equivalent damping coefficients and 

{E} is an error vector resulting from the substitution of ~c_]{Z} for 

[c 1Z}. Because the vector {Z} is random, elements of the error vector iE} 
are also random. In Eq. 4.54, there are two unknown matrices, rc] and {E}. 
[It is therefore possible to choose re] in such a way as to minimize the 
errors {E+. Because the errors are random, they will be minimized ina 
mean square sense. 


Equation 4.54 represents nm equations of the form 


n . e 
i OMNI Meet al OR ee (4.55) 


* 
The value of CT ™ is determined from the requirement that the mean square 


error, ELE‘], be minimized. Thus, 
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2 EEC] = 0 (4.56) 


* 
Yat 





oC 


The expectation and differentiation operators in Eq. 4.56 are linear and 


thus may be interchanged. 


aE. 
ORE ey, seep (4.57) 
Y ac 


Gar 





From Eq. 4.55, the random error and its partial derivative are expressed 





as 
nm : re | 
Soe ) Cy i es (4.58) 
S=] 5S 
aE. 
contact Se (4.59) 
aC. 


Substitution of Eqs. 4.58 and 4.59 into Eq. 4.57 leaves 

m eee an 
aia ae rhe As al Ow 3 ie (4.60) 
Interchanging the order of the expectation and addition operations yields 


nm one " “9 
- “a . E[Z,.Z.] + C, E[Z] = 0 (4.61) 


Because Z,, 1S a zero-mean random variable for all r, 


Elz eZ ] = of : 
yrs LZ. 
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which 1s the covariance of the time derivatives of the normal coordinates 


for modes r and s. Furthermore, 


a ome 
EIZ.] = «§ 
Yr ZZ. 


which 1s the vartance of the time derivative of the normal coordinate for 
mode r. 


The solution of Eq. 4.67 for e . 1s 


> 


(4.62) 


* 
The selection of C. for all r using Eq. 4.62 minimizes the mean square 


3 


value of elements of the error vector E because 


PE (EC) 9@(E4) 
era | Enea 
Gre aCe 
2 £ fz 
ee ae 
= 05 5 (4.63) 
r vr 


and the variance 1s always non-negative. When the error is minimized, it 


1S approximately true that 


[C,lizt = Pe iz (4.64) 
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Substitution of Eq. 4.64 into Eq. 4.52 completes the uncoupling of the 


equations of motion. 
oo P a * A o ? . 7 
miZe + PO Z)o> me oe ea talses (4.65) 
Equation 4.65 represents nm independent differential equations of the form 


2 eee 
Mee tO een 


Z 
oni iesteme er 


a col) 4 Tp} (4.66) 


4.4.2 Solution of the Uncoupled Equations of Motion 


The solution to Eq. 4.66 is obtained by using the convolution 


integral (also referred to as Duhamel's integral) to obtain an expression 


th 


for the random normal coordinate of the r~ mode. First each term of Eq. 


4.66 1s divided by m.,: 


* 


2 
Z 
e et = 


co") 3 Tepy (4.67) 


ee fale = 
Zh 
Y m Y 


= he 


Next, two new variables are introduced. Let Cn be the damping ratio for 


mode r, defined by 


* 
and let P(t) be the random modal force for mode r defined by 


eel. co py mee Tee S ale) 
i i nF bey onk Prk’) (4.69) 


The substitution of Eqs. 4.68 and 4.69 into Eq. 4.67 leads to 
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Z, + 26,02, + wil, = P(t) (4.70) 


Equation 4.70 is of a form commonly encountered in structural dynamics pro- 
blems. The solution for Z(t) may be written using the convolution inte- 


gral. 32 
ic 


Z(t) = f PRC nel tes) ogy (4.71) 


a 


The function, h(t), which 1s commonly known as the impulse response func- 


tion, is defined by 


h(t) =< 0 : ta. 0 


HE WT 
ge et sin (Vl - CWT) ; ee 
Ol shCUEL 5 


w Vi - 22 
r Yr 


The variable = 1s simply a dummy time variable used in the integration pro- 
cess. The lower limit of integration, -~, implies that the random process 
has been going on so long that all transient starting effects have been 


damped out. 
4.5 Derivation of Response Spectral Density Functions 


At this point, it is worthwhile to back away from the details of 
the problem and look again at the overall picture. The goal is a mathe- 


matical model relating the riser response spectra, Sy y (2). S. : (2), and 
ni nj 00 

Sy My (2), to the sea surface elevation spectrum, S_(%) So far, a rela- 
ni nj 


tionship between the sea surface elevation spectrum and the water particle 


kinematic spectra has been developed with the use of linear wave theory. 
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This link, which is given by Eqs. 3.48, 3.54, and 3.56, is valid for both 
static and dynamic problems. 

In this section, expressions relating the kinematic spectra to 
the response spectra are derived in four steps. Equation 4.42 gives the 
random force intensities, {P(t)}, in terms of the random fluid kinematics, 
{U(t)} and {U(t)}. The random modal forces, (P'(t)}, are expressed as 
functions of the random force intensities by Eq. 4.69. The random norma] 
coordinates, {Z(t)}, are determined from the modal forces by Eq. 4.71. 

And the random node deflections, {Y(t)}, may be calculated from the normal 
coordinates by means of Eq. 4.47. With these four relationships, it is 
possible to develop formulas relating the corresponding spectral densities 


in a manner similar to that used in Chapter 3. 
4.5.1 Spectral Density of Random Wave Force Intensities 


The expression for random wave force intensity for the linearized 
dynamic problem, Eq. 4.42, is very similar to the equation for random wave 


force intensity for the linearized static problem, Eq. 3.72. 


Pix tf) = CAB cail(xt) + C_U(x,t) (4,42) 
POGti = 67/8 onbexst) + CU(x,t) (3272) 


The only difference between the two force equations is that the drag term 
of the static equation is a function of the standard deviation of the ab- 
solute water particle velocity, ote while the drag term of the dynamic 
equation 1s a function of the standard deviation of the water particle velo- 


city relative to the structure, Oy The force intensity spectral density 





oF 


function for the dynamic problem could be derived in the same fashion as was 
used to derive Eq. 3.78 for the static problem. However, an easier derivation 
is possible by using Eq. 3.78 as a guide and replacing O14 by Oy wherever it 
occurs in that equation. When this is done, the spectral density of the 


random wave force intensity for the dynamic problem is given by 


28 8 
(ay = es 2c. oc 4 190 ,¢\/8 (oj -o, ) 
Pp Pie U T ure ae ua T Nee Vee 


nk 


+ wee | Se (2) (4.73) 


Unk ene 


4.5.2 Spectral Density of Random Modal Force 


th 


The random modal force for the r”~ mode is given by Eq. 4.69. By 


means of Eq. 4.69, the correlation function of the modal forces for modes 


r and s may be written 


* * 
Retpt(o) = ELPL(t)Pe(t + 7) 


i 5 oft)p (ty Fg 
nk= k=] 


_y?nk Pnk'®) ae nett + t)] (4.74) 
en 


ne 


In Eq. 4.74, the summations may be written as a double sum and the order 


of summation and expectation may be interchanged to give 


epee 


. (r) (5) 
; L one o ne El , bt ) Pett a t) ] (a7 5) 
Y Ss n nc= 


=] 


i 
Mm M\— 


But the expectation term in Eq. 4.75 is simply the correlation function of 


the random wave force intensities at nodes nk and né. 
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Therefore, 
(1) (4.76) 
Pak ng 


The spectral density function of the modal forces is the Fourier 


transform of the correlation function, Eq. 4.76, 


1 1 rf 
2) il y  } oO %6 
oe Py 2 J me nk=1 né=] K ne 


(4.77) 


Pak’ ng 


In Eq. 4.77, the order of integration and summation may be interchanged 


: 1 Ff 4 Nos (s) 
S * *(9) = panes , > 
eA me nk=1 né=1 Te 


] { -12% 
oa R (r) e dr] (2573) 
on Me Pat 


The term in brackets in Eq. 4.78 is the spectral density function of the 
random wave force intensities, . p (2). Therefore, the one-sided spectral 
nk ng 


density of the modal forces is 


sptpnta) = Te bolt) ol) sy oe) 
rs Ma nk=1 né=1 nk n& 
4.5.3 Spectral Density of Random Normal Coordinates 
th 


Equation 4.71 expresses the random normal coordinate of the r 
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classical normal mode as a function of the random modal force for that mode. 
In this section an expression relating the corresponding spectral densities 
mseaerived. 

In the usual fashion, the correlation function of the normal co- 
ordinates for modes r and s is written 

fae = E[Z.(t) Z(t + t)] (4.80) 
feom Eq. 4.71, Bie + t) is given by 


cet 
* 
Z(t +r) =f Pele helt + «= &5) de, (4.81) 


where Eos like Ey: 1S a dummy time variable used in the integration process. 
When Eqs. 4.71 and 4.81 are substituted into Eq. 4.80, there results an ex- 


panded formula for the correlation function. 


iL 


Re) > et J Prle)ny ayes 
Vere 
f P(Ep)h.(t a Dieee Dy de (4.82) 


-_ X& 


The spectral density function of the random normal coordinates is 


the Fourier transform of the correlation function. 
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Sua = Ls tf h(t - &))dé, 


r s 
T==_0 1a 
Ctr 
f p- g(Eo)h, (ttt- E,)dé, le int y (anoor 
ane 


Evaluation of Eq. 4.83 is simplified by a convenient transformation of 


coordinates. Let 


E4 Ee Ey (4.84a) 


Eq = 1h ay 69 (4.84b) 
where b3 and Eqs like Ey and Eos are simply dummy time variables. The use of 


Eqs. 4.84 transforms Eq. 4.83 to 


- -{2QT 
(- f Pet + x = eg dig (eq)deg) Te” ta (4.85) 
7 aes 
The terms in Eq. 4.85 may be more conveniently arranged by moving the expec- 
tation operator inside of the second and third integrals and reversing the 


integration limits of those integrals. 
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© °°) foe] 


ak Lf J o%-or 
sale ok [Prt - c,)Pe(t +1 - &,)] 
rs ae so, ee 

TH-0 E3°0 E,=0 


nn(Eg)he(Eg)e  “degdegdr (4.86) 


In Eq. 4.86, the expectation term is a correlation function, 
* * _ ae 
ELPH(t - Eg)Pg(t + 1 = Eq)] = RptpA(e = ky + E3) (4.87) 
and the exponential term may be decomposed into three components. 
ES Ete 5) ee ye len (4.88) 


The introduction of Eqs. 4.87 and 4.88 into Eq. 4.86 results in 


8 
8 


co 


~ . ie _ 
7) a, { if | Epep oe 3 aaeea! 
{T==-0 E3=0 E470 


-i2(t + & - &,) ‘ 
e ee) eee el 


[ho (Ege “4 ]drde dé, (4.89) 


Because h(t) = 0 for t < 0, the lower limits of integration with respect to 
b3 and Eq may be changed to -~ without changing the value of Eq. 4.88. When 
this is done, the two integrals with respect to 54 and Eq may be recognized 
as the Fourier transforms of the impulse response functions. From the theory 


of random vibrations,*” the Fourier transform of h.(&,) is 
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f ngleqieT 84 dey = H,(0) (4.90) 


where 


H.(2) = = (4.91) 
: wf fare 12¢w.0 


The function, Ho (2), given by Eq. 4.91, is commonly referred to as the 


frequency response function or transfer function of the system. The integral 


with respect to £3 Ts 


ioe) 


I h.(Egle°3 dé, = A(ay (4.92) 


Yr 


where the bar in H (2) denotes the complex conjugate of H (2). Substitution 


of Eqs. 4.90 and 4.92 into Eq. 4.89 leaves 


ae ee *( ge, 
Sz z (2) = a Jom Rp Ge ana Eq)e dt 
nes 
: H (2) H. (2) (4.93) 


The integral in Eq. 4.93 is the Fourier transform of a correlation function, 


which, by definition, is a spectral density function. Therefore, 
C See 
i Spipe (a) H (2) H. (2) (4.94) 


and the one-sided spectral density function of the normal coordinates is 


given by 


Sy z (a) = Sp*pr (2) H (a) H.(2) (4.95) 


rs 
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4.5.4 Spectral Density of Response 
Spectral Density of Node Deflections 


Equation 4.47 expresses the random node deflections, {Y(t)}, in 
terms of the random normal coordinates, {Z(t)}, and the eigenvector matrix, 


[¢]. For a single random node deflection, this relationship is 


nm 
LG) 4 pt) Z(t) (4.96) 


Equation 4.96 states that the random deflection of node ni is the sum over 
all modes of the products of the node ni deflection in each mode shape 
multiplied by the random normal coordinate of that mode. In effect, the 
normal coordinates are weighting factors which assign the relative importance 
of the natural modes of vibration in the riser's response. 

Equation 4.96 is in a form convenient to use in deriving an ex- 
pression for the spectral density of the node deflections. As usual, this 
derivation starts with the correlation function of the random deflections at 


nodes ni and nj. 


Ry y (7) = ELY«(t)¥ .(t + r)] (4.97) 


ni nj nJ 
Equation 4.97 may be written in terms of the normal coordinates by substi- 


tuting Eq. 4.96. 


(r)7 (ty FalS) z(t + 2) (4.98) 
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In Eq. 4.98, the order of summation and expectation may be interchanged, 
the two sums nested, and all deterministic variables removed from the expec- 
tation operation. When this is done, the result is 

nm 


nm 
it) = (r) ,(s) 
Ying ; a of oni tng FLZ.At)Z,{t + 7) (4.99) 


Substitution of Eq. 4.80 into Eq. 4.99 produces 


nm nm 
C ie 

Ry oy tr) =) ane? 
ni nj r=] s=] 


(s) (4.100) 
One ok (1) 
ni ‘nj LZ. 

The spectral density function of the random deflections at node ni 
and node nj is the Fourier transform of the correlation function, Eq. 4.100. 


. 1 rf nm nm (rs) -72t 
Sy oy (a) = ef SF Lo elsde 5 (2) ota (4.101) 
ake i ni 12 


The order of tntegration and summation in Eq. 4.101 may be interchanged to 
give 


co 


nn nm é 
= ] ~7LT 
§ GQ =) eee ee { Te Ge he (4.102) 
ant os gay (UL Tis Qn LZ. 


=O 


Finally, the integral term in Eq. 4.102, being the Fourier transform of the 
correlation function, is the spectral density function of the random normal 
coordinates r and s. Therefore, the one-sided spectral density function 


for the node deflections is 
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nm nm 
S (2) a eee) Sia) (4.103) 
Yaa ng releset,) eae! Ze, 
Equations 3.48, 4.95, 4.73, 4.79, and 4.103 provide a mathematical link 
between the input surface elevation spectral density function and the out- 


put node deflection spectral density function. These five equations may be 


combined into a single expression. 


eae . git T a ol), H(S) H 
se net ° nk ne t me “Vik vne 
ae = le, aes ) + 0%] 
Ink (®) Ine@)} Syn (2) (4.104) 


The power transfer function of the mathematical model is the entire term 
within the braces in Eq. 4.104. 

Before deriving the bottom rotation and bending moment spectral 
density functions, it is worthwhile to consider the nature of the solution 
indicated by Eq. 4.104. There are two types of variables that make the 
solution a trial and error process, because they cannot be immediately 
determined. One of these variables is the standard deviation of the rela- 
tive water particle velocity, A at the nodes. Because the relative 


water particle velocity is a function of the unknown random structural 


velocity, Y ae values must be assumed for [~o oy and updated on each iteration 
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until convergence occurs. The second variable which changes in each itera- 
tion is the damping ratio, Ces which is used to calculate the frequency 
response function, H (2). Equations 4.68 and 4.62 show that Ch is a func- 
tion of the covariances of the first time derivatives of the random normal 
coordinates, of 5 which are also unknowns. Equations 4.68, 4.62, 4.53, 
and 4.25 show fe = also depends upon [-oy J. 

Thus, in the calculation of the random response, it is necessary 
to assume values for [oy J and [0$5] and use a trial and error process in 
which [ oy J and [o$5] are calculated in each cycle of iteration. The itera- 
tions are continued until the calculated values of the two matrices agree 


satisfactorily with the assumed values. 
Spectral Density of Bottom Rotation 


For the dynamic problem, the bottom rotation bears the same rela- 
tionship to the deflection of the first node as it does for the static pro- 
blem. Therefore, Eq. 3.93, which was derived for the static problem, is 


also applicable to the dynamic problem. 
Spectral Density of Random Dynamic Bending Moments 


In Chapter 3, the following expression was derived relating the 


random bending moment at node ni, M .(t), to the random node displacements. 


ni 


| ni+] | 
Matt)? = oe Jeomeni (3.98) 
where [J] is a transformation matrix defined by Eqs. 3.97. As indicated by 


Eq. 4.96, the random node deflections are a function of the classical normal 
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mode shapes and random normal coordinates. The substitution of Eq. 4.96 


into Eq. 3.98 yields 


ne nm (r) 
Met) ee Deen elt (4.105) 
ni Arai! nisnk oy nk “r 
or, in matrix form, 
OL) ae 8 ree col epee cs AL Re (4.106) 


Let [B] = [J][#] be an n by nm transformation matrix relating random nodal 


moments to the random normal coordinates. Then 


{M(t)} = [B]{Z(t)} (4.107) 


or for a given nodal moment, 
nm 


Ha a) 8 


= Ls Bai eeelt) (4.108) 


naar’ 

Equation 4.108 may be used to write the correlation function of 
the random moment at node ni with the random moment at node nj. 

R TR co (ee Gee ai] 

Ming ni nj 

m nm 

ob Bajsstst 7) (4.109) 
In Eq. 4.109, the summation operations may be combined and the expectation 


operation interchanged with the two summations because all three operations 


are linear. This leads to 
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nm nm 


(tea ine OTB clit. oR (=) 
Maing r=] seq Mer nds ZZ. 


(4.110) 


The spectral density function of the moments at nodes ni and nj 


is the Fourier transform of the correlation function. 


Lf ef (rye Md (4.111) 
5 (2) = 1 f Beebe tje dt 4.111 
Maing dn } y=] $21 nisr nj,s ZZ, 


Interchanging the order of summation and integration in Eq. 4.111 yields 


nm nm 


] =-7 Qt 
) B. BB. . [5 | R (t)e dt] (Aeizy 
Maing eT een hie (ijieSaeen ZZ. 


NW 


S 
The expression in brackets in Eq. 4.112 is, by definition, the spectral 
density function of the random normal coordinates for modes r ands. It 


follows that the one-sided spectral density function for bending moments is 


nm nm 


(Q) =~). ) Bee e: 
Mang pelea Ee SINSs 


S (4.113) 


ze 
Equation 4.113, together with Eqs. 4.95, 4.79, 4.73, and 3.48, 
forms a mathematical chain of relationships linking the input random sea 
surface elevation spectrum to the random nodal moment spectra. From the 
Similarity between Eqs. 4.113 and 4.103, it follows that a single expression 
for bending moment spectral density as a function of the sea surface eleva- 
tion spectral density can be expressed by Eq. 4.104 with See substituted 


Tor pt) and Bj : substituted for je 
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4.6 Damping 


Although it may not be obvious, the response spectral density 
functions of the riser depend upon the damping matrix, [c]. That this is so 
may be seen by considering Eqs. 4.104, 4.91, 4.68, 4.62, 4.53, and 4.25. 
These equations show that the response is a function of the frequency re- 
sponse functions, H (2), which depend upon the damping ratios, Cpe which 
in turn are functions of the diagonalized modal damping coefficients, ce 
The diagonalized modal damping coefficients are calculated from the modal 
damping matrix, [CJ]. which is a function of the modified damping matrix, 
[c]. Because the modified damping matrix is a function of the damping 
matrix, [c], the response spectral density function is also a function of 
the damping matrix. 

The damping matrix was introduced in the damping force term of 
Eq. 4.4 without any indication as to how the coefficients which constitute 
the matrix may be determined. In this section, a rational method of formu- 
lating the damping matrix is developed. 


Consider the modal damping matrix as defined by Eq. 4.53. 
[c.] = [e]'[elfo] (4.53) 


The modal damping matrix results from a coordinate transformation which is 
used to uncouple the simultaneous differential equations of motion. If 
[c] is a diagonal matrix, that is, if the normal modes are orthogonal with 
respect to the modified damping matrix, [c], then the equations of motion 
are indeed uncoupled by the transformation. The modified damping matrix 


is given by 
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Cc] = [cl + cB Poy J (4.25) 


Substitution of Eq. 4.25 into Eq. 4.53 leads to 
3 I . 8 lee 
[C.J = (Ce) [elle] + cA/* [o]' boy JLo] (4.114) 


The second term on the right-hand side of Eq. 4.114 is not, 1n general, a 
diagonal matrix, because there 1s no reason for the normal modes to be 
orthogonal with respect to the matrix of standard deviations of the rela- 
tive velocities. Therefore, the modal matrix will not, in general, be a 
diagonal matrix, and the diagonalization procedure described in Section 4.4.1 
must be used to uncouple the equations. The error resulting from this 
diagonalization may be reduced somewhat, however, if the modal damping 
matrix is as nearly diagonal as possible, that is, if the diagonal elements 
of Koel are large with respect to the off-diagonal elements. The diagonal 
elements of the modal damping matrix can be strengthened with respect to the 
off-diagonal elements if the first term on the right-hand side of Eq 4 114 
is a diagonal matrix. 

It may be concluded that a logical criterion governing the formu- 
lation of the damping matrix, [c], 1s that the normal modes be orthogonal 
with respect to the damping matrix. In section 4.6.1, the conditions under 
which this occurs are explored. 

4.6.1 Conditions Under Which Normal Modes Are Orthogonal With Respect 
to the Damping Matrix 

The conditions necessary for the normal modes to be orthogonal 


with respect to the damping matrix have been demonstrated by Caughey?; and 


Ma 


are repeated here without proof. Caughey shows that simultaneous equations 


of motion expressed by 
(-IJip} + [A]{p} + (B]i{p} = {F(t)} (4.115) 


are uncoupled by the same transformation which diagonalizes the undamped 
system if 
» N-1 ; 
iL 
is =e “al (4.116) 
j=l 2=0 

Equation 4.116 is a sufficient condition, but not a necessary condition. 
In Eq. 4.116, N is the order of the matrix [B], the subscript 1 on the matrix 
[B] means that only one root of [B] is taken, and a0 is a linear coefficient. 

Caughey's findings may be applied to the riser problem, for which 


the linearized equation of motion is Eq. 4.31, 
[m,~]¥} + Cclt¥} + [K](¥} = {P(t)} (4.31) 


where {P(t)} is defined by Eq. 4.42. A form similar to that of Eq. 4.115 
is obtained by premultiplying Eq. 4.31 by areal When this is done and 
Eq. 4.25 is substituted for [c], the result is 


11} = fe ene + aye Pied Psi 


+ EmIo'eKI Y= Om 'eP(t)) (4.117) 


From Eqs. 4.116 and 4.117, it follows that [oJ fclfe] is a diagonal matrix 
rT 


lz 


fe o} N-1 ge 
-] gal 
melt Lelie a8] (4.118) 
: iO op ae 


where 


[B] = Cm, J'EK] 


Premultiplying Eq. 4.118 by tm~] yields the desired expression for the 
damping matrix. 
co N-] : 
¥ £ | 
fee = Emel). cele: (4.119) 
j=1 2-0 


For illustration purposes, consider the case where the coefficient 
aie 1S non-zero only when j is 1 and £ 1S zero or one. From Eq. 4.119, the 


damping matrix for this case is given by 


Lc] [~m~] CHA e + a Be 


ayot-m,~l + a4, LK] (4,120) 


Equation 4.120 describes Rayleigh damping, which 1s frequently assumed in 
dynamic analysis in order to insure that the orthogonality transformation 
will uncouple the equations of motion. Because the index j in Eq. 4.119 has 
no finite upper limit, Eq. 4.120 is only one of an infinite number of ways 
in which the damping matrix may be formulated such that the normal modes 


are orthogonal with respect to it. 
4.6.2 Formulation of the Damping Matrix for the Marine Riser 


Consider the schematic representation of a marine riser shown in 
Fig. 4.6, in which damping is represented by a series of dashpots. The 


damping forces may be divided into two classifications--environmental damping, 
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which depends upon the velocity of the structure with respect to the sur- 
rounding medium, and internal friction damping, which depends upon the rela- 
tive velocity of one part of the structure with respect to other parts of 
the structure. In Fig. 4.6, environmental damping is represented by the 
dashpots labeled Coe and internal friction damping is indicated by the 
dashpots labeled Ces which connect adjacent nodes. Although they have been 
omitted for the sake of clarity, dashpots representing internal friction 
damping could be shown connecting each node to every other node. 

When the equations of motion of the riser shown in Fig. 4.6 are 
written in matrix form, the damping matrix is composed of elements which 
are combinations of the coefficients Cy and Ce- The environmental damping 
coefficients appear only in the diagonal elements of the damping matrix, 
but internal friction damping coefficients contribute to all elements of the 
matrix. Thus, environmental damping is similar to the first term of Eq. 
4,120 in the sense that both contribute only to the diagonal elements. In 
like manner, internal friction damping is similar to the second term of 
Eq. 4.120 in the sense that both contribute.to all elements of the damping 
matrix. 

Let 8. be the damping ratio for mode r defined by 


= eal (rata eclr) 7 
aoe Ornate ton cio a (4.121) 


and suppose that [c] is formulated according to Eq. 4.120. Then 


aicuee Te (ay toa ema y + aco eK 3) (4.122) 
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Utilization of Eqs. 4.39 and 4.46 in Eq. 4.122 leads to 


] | 2 
8B = gs— (a,,m_ + ajywm) = i 
r 2mw, 10-e llre eu), Z 








(4.123) 


Equation 4.123 shows that if the damping matrix 1s proportional to the mass 
matrix alone, then the damping ratio decreases with increasing natural fre- 
quency, and higher modes are damped less than lower modes are. Further, 

if the damping matrix is proportional to the stiffness matrix alone, then 
the damping ratio increases as frequency increases, and higher modes are 
damped more than lower modes are. 

It can be shown that the damping matrix, [c], can be completely 
specified by selecting the modal damping ratios. Let [~g~] be an nm by nm 
diagonal matrix of modal damping ratios and let fw] be an nm by nm dia- 
gonal matrix of natural frequencies. Then, if [c] is chosen so that the 


normal modes are orthogonal with respect to it, [8 _] is given by 
] at = Zor = 
pe uel lcliely = Bel bey ae) (4.124) 
e 


Premultiplication of both sides of Eq. 4.124 by [o], postmultiplication by 
el, and utilization of Eq. 4.45 leads to 


[c] = amlelbsJb«tte! (4.125) 


Equation 4.125 reduces the problem of selecting n@ damping coefficients to 
one of selecting nm damping ratios. 
Consider now the selection of the damping ratios, Bis Y = Veo 


. «5 nm. The damping ratios may increase with mode number, they may 
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decrease with mode number, and they may be constant for all modes. 

The case in which the damping ratios increase with mode number 
occurs if damping is proportional to the stiffness matrix alone, a condi- 
tion which has been compared to internal friction damping. If the last 
term of Eq. 4.120 1s used to represent internal friction damping, and if 
internal friction damping predominates over environmental damping, then 
the damping ratio for higher modes will be greater than it is for lower 
modes. 

The case in which damping ratios decrease as mode numbers in- 
crease occurs 1f damping is proportional to the mass matrix alone, a situa- 
tion which has been compared to environmental damping. If the first term 
of Eq. 4.120 is used to represent environmental damping forces and if 
environmental damping predominates over internal friction damping, then 
higher modes will be damped less than lower modes are. 

A third alternative is for the damping ratio to be constant for 
all modes. Because Eq. 4.120 contains only two coefficients which may be 
varied, a constant damping ratio may be specified for only two modes if 
Rayleigh damping 1s assumed. However, a constant damping ratio for all 
modes may be specified by the more general combination of mass and stift- 
ness matrices indicated by Eq. 4.119. A more direct method of specifying 
constant damping ratios is by substitution in Eq. 4.123. 

Based on the foregoing discussion, an argument for the specifi- 
cation of the damping ratios may be made. For the marine riser moving 
through the water, the environmental damping force is simply the drag 
force of the water which resists the motion of the riser. In the formula- 


tion of the equation of motion for the riser, Eq. 4.8, this drag force 15s 
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accounted for by considering wave structure interaction. The drag force 
caused by wave motion, which is an exciting force, and the drag force caused 
by the riser's motion through the water, which is a damping force, are both 
included in the single term C-{V| VI}. Because the environmental damping 

is treated as an applied force in Eq. 4.8, the damping term, [c]{¥}, repre- 
sents only the resistance to riser motion caused by internal friction. 

Both internal friction damping and environmental damping are included in 
the modified damping matrix, [c], which results from the linearization of 
the equation of motion. This is shown best by Eq. 4.25, the first term 

of which, [c], represents internal friction damping and the last term of 
which is a modified drag coefficient. 

The damping ratios, [~8 ], which determine the damping matrix, 
[c], may be considered as internal friction damping ratios. Given the 
argument that internal friction damping is similar to that component of 
Rayleigh damping which is proportional to the stiffness matrix, it appears 
unlikely that 8. decreases with mode number. On the contrary, it 1s 
logical to assume that B. increases with mode number. But if the method of 
analysis permits B. to increase too much with r, significant contributions 
of higher modes may be damped into insignificance. It is reasonable, there- 
fore, to select [c] such that g. 1S constant for all modes, with the expec- 
tation that such a procedure will make any errors in higher mode contribu- 
tions tend to be conservative ones. 

There are precendents in the literature for using a constant 
damping ratio for all modes. Caughey?? gives an example of this type of 
damping and refers to it as linear hysteresis damping. Hart°* uses a damping 
ratio of .02 for all modes in considering the response of steel building 


frames. 


Ly 


In this thesis, the damping matrix, [c], will be determined in 
such a way that B is constant for all modes. When this assumption 1s made, 


Eq. 4.124 simplifies to 


[o]'[c]lo] = 2m,ebw.] (4.124a) 


and Eq. 4.114 becomes 


[C,] = 2mefw] + ca/® tel" C-os Ite] (4.114a) 


The first term of Eq. 4.114a contributes only to the diagonal elements of 
[CJ and represents modal damping due to internal friction in the riser. 
The second term represents damping due to fluid drag forces resulting from 
the relative velocity between the riser and the surrounding water. 
4.7 Variances and Covariances of Time Derivatives of Random Normal 
Coordinates 

The simultaneous equations of motion are uncoupled by using the 
orthogonality of the normal modes with respect to the mass matrix and by 
employing a diagonalization technique to eliminate the off-diagonal ele- 
ments of the modal damping matrix, (CJ. As indicated by Eq. 4.62, ele- 
ments of the diagonalized modal damping matrix are functions of the vari- 
ances, of and covariances, of # Ot The First Lime denivativesvore che 
random normal coordinates. The variance or covariance may be calculated by 


integrating the corresponding spectral density with respect to frequency. 


2 
5 = ji S54 (a) da (4.126) 
0 CS 
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In a fashion similar to that used in Section 3.3.4 to derive Eq. 3.56b, 


it can be shown that 
S5 5 (Oe ea G Sz 7 (2) (4.127) 
rs r 
Combining Eqs. 4.126 and 4.127 leads to 


ae 2 
oF 5 = fo Sz 7 (2)d2 (4.128) 


Equation 4.128 may be used to calculate the variances and covariances used 


to diagonalize the modal damping matrix. 
4.8 Calculation of Standard Deviations of Relative Veiocities 


In computations for the random riser response, the standard devia- 
tions of the relative water velocities with respect to the structure, [~oy I, 
are needed in two instances. These standard deviations are first used to 
calculate the modified damping matrix, [c], according to Eq. 4.25. In this 
instance, the relative velocity standard deviations are a measure of the 
environmental damping of the dynamic structural response by the water. 
Secondly, the relative velocity standard deviations are used in Eq. 4.73 
to calculate the spectral density functions of the linearized force inten- 
Sities at the nodes. 


The random relative velocity at node nk, as defined by Eq. 4.5a, 


Vist) = Ut) - ¥ tt) (4.129) 
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The autocorrelation function of the relative velocity 1s 


yo gt = SEUNG EIN A Geer I me yee Yale) 


(ine ae ele aay (4.130) 


ak 


Expansion of Eq. 4.130 and interchanging the order of addition and expec- 


tation operators leads to 


a ELU, _(t)u nkbt * td] - SPT 0 lio 2 7) 
nk nk 


: EU. (t¥ p(t + t)] + Bt teyce + +)] (4.131) 


Each of the four terms in Eq. 4.131 is an autocorrelation or cross correla- 


tion function. Therefore, 


R° ° R° ° ( = Re e ¢ R° ° (a 
vaenk Cane Meena valent 
+ Re tees) (4132) 
Japan 


The spectral density of the random relative velocity at node nk is obtained 
by taking the Fourier transform of Eq. 4.132. Because the Fourier trans- 
form of each of the four terms on the right-hand side of Eq. 4.132 is a 
spectral density function, there results 


SH (0) Se oe eee Soren nt 0), = hoe ecm) 
Yak’ nk Yak Unk Ulf Tak unk 


yy (22) (Aeleo) 
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Finally, the variance of the relative velocity at node nk is the integral 


of the spectral density with respect to frequency, 


of - [ Sy J (ayaa (4.134) 
0 n n 


and the standard deviation is the square root of the variance. 
ety = af (4.135) 
V V : 


Consider the four terms of Eq. 4.133 which contribute to the spec- 


tral density of the relative velocity. The first term, Si U (2), is the 
nk nk 
spectral density of the random water velocity and is given by Eq. 3.48. The 


last term is readily calculated because 


(a) = as, y (a) (4.136) 


nk nk 


Ss e 

eink 

where Sy y (2), the spectral density of the random deflection at node nk, 
nk nk 

1s determined from either Eq. 4.103 or 4.104. The two cross spectral den- 


sity terms remain to be determined. By analogy to Eqs. 3.54 and 3.56a, it 
follows that 


Soo A(R) = Ses = (2) (4787) 
Orel teal 

and 
Se. oe (= aS oe oem) (4.138) 
ane Ynkenk 


A general property of cross spectral density functions is that they are 


hermitian,2’ that is 
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Sy ( (2) = Si) y (2) (4.139) 
nk “nk nk nk 


where the line over the spectral density function denotes the complex con- 
jugate. Therefore, 
Sy ( (2) = -i2 Sf y (2) (4.140) 
nk “nk nk nk 
and it is only necessary to derive an expression for Si y (2) to determine 


nk nk 
the two cross spectral density functions of Eq. 4.133. 


4.8.1 Cross Spectral Density of Water Velocity and Node Displacement 


The random displacement of node nk may be written in expanded form 
by using Eqs. 4.96, 4.71, 4.69, and 4.42. 


e 
1 ate ronal eG) 8 . 
Y(t) = — | 6 Pe, Cay ona) 
nk Mm, a nk nee nz uV 1 Vie ne] 
] 


+ CU p(E,)) Ih, (t-e, de, (4.141) 


The cross correlation function of the water velocity and deflection at 


node nk is 


Go lh) ACE Se A] (4.142) 


nk nk 


Equation 4.142 may be expanded by the introduction of Eq. 4.141 for the node 


displacement. 
nm CNT 
: Z i ] mee) (ie) 8 
R (r) = ELU (t) — } 6 f yoy? (Ca [= 
Biel a nk ls Sen lk ip es ie J 
lias 


mH onesey + CU p(E,) h(t + t- &))de,] (4.143) 
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In Eq. 4.143, because the expectation operation, summations, and integra- 
tion are all linear operations, their order may be interchanged. When 


this is done and the substitution, Eo = tet ate Ey 1s made, there results 


— 1 © tr) Ff ,fr) 2 
2 an 1 *nk wayne tC TV, 


nk e r= 
[ebb (tvb,plt + + ~ 65)IR(6y)d6, 
0 


Paes f JU elt = aco 
0 


h (E> )dé, } (4.144) 


By definition, the impulse response function, h.(E5)» is zero for Eo < 0. 
Therefore, the lower limits of the integration in Eq. 4.144 may be changed 

to -~ without changing the value of the equation. The two expectation terms 
in Eq. 4.144 are recognizable as the correlation functions, Ri U (t-&5) and 


Ry «| (t-&) respectively. If &, = 1t - &, Eq. 4.144 may be written 
Une 2 3 2 


e (r) = 1 g(t) Y ot) rca/® oj | R- . (é ) 
Unk nk Me pz? MK nteq ie INS he _ Unkene 3 


t°.¢) 


RelEp)86> + Cy f Ry ij C6 Bp(Ep)¥ey] (4.145) 


The cross spectral density function of the water velocity at node 
nk and the structural displacement at node nk is the Fourier transform of 


the cross correlation function. 
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See = - { ye 18tg, (4.146) 
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In Eq. 4.146, the terme may be written AS. | e7 1883 by the de- 


finition of E3° Substitution of Eq. 4.145 into Eq. 4.146 produces 


~ 


lau (r) © o(r) 
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oe ak m, en Es y nk ya] % 
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aes 


co) in, 
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“128, 
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Integration of Eq. 4.147 with respect tog, yields 
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Finally, integration of Eq. 4.148 and substitution of Eq. 3.54 results in 


the following expression for the one-sided spectral density. 


1 a (r) f  .(r) 
Sy le) — J) H (a) ¢ 6 
Unk nk Mm r=] ° nk nl=1 nt 
8 
ccaE oj Pee Le) (4.149) 
u\/ 7 Vie a Ui ene 


4.8.2 Spectral Density Function of Relative Water Velocity 


The spectral density function of the relative water velocity may 
be written by substituting Eqs. 4.136, 4.137, and 4.140 into Eq. 4.133. 


Sy y (oie = ss" 


So. 16S: (2) 
nk “nk Dak enk Unk! 


nk nk 


+ 72 Si y (2) + a’sy y (2) (4.150) 
nk nk nk nk 


The second and third terms on the right-hand side of Eq. 4.150 are complex 


conjugates, the sum of which is simply twice the real part of either term. 


“19S y (2) + i2 Sty sy (a) = 2 Re[-i2Si, y @)] (4.151) 
nk“ nk nk nk nk nk 


The real part of these terms may be evaluated from Eq. 4.149. 


nm n 
Re[-i2S- (2)] = a alt) : gt) {Re[H (a) ] 
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In Eqs. 4.151 and 4.152, Re indicates "real part of," and Im indicates 
"imaginary part of." By using Eqs. 4.151 and 4.152, Eqs. 4.150 may be 


written in the following form. 


2 
Set) = | ohne (co) eas (2) 
vat ak Uk unk Yak’ nk 


nm n 
ae) (0) ) p(t) {C.2- Re[H(2)] 
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. cal oi aC ees aac) (4.153) 
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Because each of the variables in Eq. 4.153 has been previously 
derived, Eqs. 4.153, 4.134, and 4.135 may be used to calculate the standard 
deviation of the relative velocity. The appearance of this standard devia- 
tion on the right-hand side of Eq. 4.153 necessitates the use of an itera- 


tive procedure in response calculations. 
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Chapter 5 


THE DYNAMIC RESPONSE OF A MARINE RISER TO COMBINATIONS OF RANDOM 
WAVE FORCES, RANDOM TOP OFFSET, AND DETERMINISTIC CURRENT FORCES 


5.1 General 


So far, the only forces which have been included in the riser anal- 
ysis are those caused by random waves. It is not uncommon that a marine riser 
installation be subjected to a steady current force as well as random wave 
forces. Furthermore, unless the surface support platform is fixed, there is 
almost always some offset of the riser top, which results from lateral mo- 
tion of the surface support platform. The introduction of current and top 
offset produces additional forces on the riser which, when combined with the 
random wave forces, lead to a response which is different from that caused 
by random waves alone. In this chapter, a method is developed for determining 


the total response to random waves, steady current, and random top offset. 
5.2 Response to Random Wave Forces and Steady Current Forces 


5.2.1 Equation of Motion 


Consider a marine riser which is subjected to random wave forces 
and a deterministic steady current which varies in some known manner along 
the riser axis. For the present discussion the top offset is assumed to be 
zero. Let U. (x) be the steady current velocity and {uo} be an n by 1 vec- 
tor of current velocities at the riser model nodes. The total relative wa- 


ter velocity due to current, waves, and structural response is given by 


T27 


V(xt) = a (x) + Ut) =F (xt) 


S Un (x) + WV (x,t) (5.1) 


and the equation of motion of the riser model (Eq. 4.8) becomes 


[-m.] iv} + [ec] (¥} + [kK] (¥} 


= C. AE C. {V} (5.2) 
While this equation of motion is somewhat similar to Eq. 4.8, there 
is one highly significant difference. The velocity parameter, V (Kot) asa 
random variable whose mean is nonzero. In fact the mean value of the random 


velocity is exactly equal to the current velocity, Un (x). 
© I Gao) ata I (2) Mea eet 


= io (x) (553) 


Because this input parameter has a nonzero mean value, the response of the 
Structure will also have a nonzero mean value. Therefore the complete so- 
lution for the response must include a determination of the mean response 
as well as the response variance. 

Because the drag force on the riser is a nonlinear function of 
the relative water velocity with respect to the structure, the total riser 
response cannot be obtained by superimposing the separate responses to the 
Steady current acting alone and the random waves acting alone. Removal 
of the nonlinearity by the equivalent linearization technique leads to co- 
efficients which are functions of both the steady current velocity, {uo}; 


and the standard deviation of the relative water velocity due to random waves, 
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[~oy ~J. Therefore, the current and wave effects remain coupled in the solu- 
tion of the problem, and the response is expressed in terms of the mean and 
variance instead of in terms of a component due to current and a component 


caused by waves. 
5.2.2 Linearization of the Equation of Motion 


The linearization of Eq. 5.2 begins with a transformation of co- 
ordinates which is achieved by substituting Eqs. 4.9 in Eq. 5.2. The result 


is 
[-mJ{Uy + [ce] (U3 + [kK]tu} = Pm tv} 
+ ce WiVis + Eel) + Ck] (5.4) 
To both sides of Eq. 5.4, the term [c]{u.} is next added, leading to 


[-m]{0} + [c](U + ud + [K]{U} 


Tm ~iV} + C. (ivi) + [eliv) 
+ [K]{V} (oa) 


The goal of the equivalent linearization process is to replace the 
two velocity terms on the right-hand side of Eq. 5.5 with a linearized ver- 


sion of the two terms. Let 
Cc. ViVi} + [e]e¥} = {c} + [cl(V} + {E} (5.6) 


where {E} is a vector of random errors which result from the linearization, 


and {c} and [c] are matrices of coefficients which are chosen in some optimal 
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manner. Equation 5.6 is similar to Eq. 4.11 except for the vector of con- 
stants, {c}, which is introduced because of the nonzero mean nature of V (x, 2) 
By assuming that the off diagonal terms of [c] and [el are identical and 


letting [c] - [c] = [“c.~], the error vector in Eq. 5.6 may be written 


(Ey = - Co }iv} + C, (W]v|} - {ec} (5.7) 


A typical element of the error vector is 


© ® 


~ 


Chi = = "ey: Vag ts Ce Orel oa on; (5.8) 


With two unknowm coefficients in Eq. 5.8, it is possible to satisfy 
two conditions concerning the random error, En The two conditions which 


will be satisfied here are to choose Ce and C - in such a way that the mean 








nj wy 
Square error iS minimized. 
3 2 _ 
aC, E te ee = 0 (5.9a) 
nj 
ae: [Er] : 0 (5.9b) 
La 


Interchanging the order of the differentiation and expectation in Eqs. 5.9 





leads to 
9E As 
el eo | =O (5.10a) 
nj ec, 
nj 
and 
aE, 
PETE. — =] = 0 (5.10b) 
MN 9¢ 
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When Eq. 5.8 is substituted into Eqs. 5.10, the following two simultaneous 


equations in Ce and on result. 


nJ 
g rm se : =9 m 
“en E ae + Cn; E Ree = tes rg lenail! (5.11a) 
“aap EB [Saad + nj cS Ce E ag eel (Sab) 


Consider the expectation terms on the left hand sides of Eqs. 5.11. 
It has already been shown in Eq. 5.3 that the expected value of Vg is the 
Furthermore, the second moment of V4 is 


. J 
nJ 
equal to the sum of the square of the mean value and the variance as shown 


current velocity at node nj, Un 


by 


“2 _ ane 
E [V4 = E se ‘ Laue ] 


2 : : wy 
= Efu- J] + €E [2u Nes? aileron 
on Cag nj nj 
= i“ aF 04 (Z): 
nj nj 


mies introduction of Eqs. 5.3 and 5.12 into Eqs. 5.11 results in 


2 ~ ee 
CG [u cain sileeromieks = Ce Eve vie le ( ona) 
en j nj Vig nj Chj U ny ng 
and 
Cat SG eee Eee (5.13b) 
enanacns nj u nj) nj 


The evaluation of the two expectations found on the right hand 


Sides of Eqs. 5.13 is simplified by the introduction of a new random variable, 
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a, defined as 


nJ 
oe ° 
nj Cue, Niger 
nJ cory (ony 
Nee Vag 


o. e e e e .* e 
Because Vag 1S zero mean and Gaussian, it is readily seen that Vad is a 


standard Gaussian random variable with a mean of zero and a variance of 


one. From Eq. 5.14, it follows that 
+ 
Wea: Gee ee (5.15) 
nj Chj Ving nj 


e 
~ 


Another useful parameter, Onj? is defined as the ratio of the 


Current velocity to the standard deviation of the relative velocity caused 


by waves at node nj. 


(5.16) 


nj 
G .* 
From Eq. 5.15 1t 18 apparent that Ving is negative for all values of Ving 


~ x 
less than -o . and ue is positive whenever Vag is greater than -a.. 


* e Ps ox 
= + oy < 0, Vi.< =a 
nj Cj Vag nj nj nj as 
* ; 
> 0 r) Vag > ng 


Equations 5.17 are useful in evaluating the expectations on the right-hand 


meside of Eqs. 5.13. 
By the 


Consider the term on the right-hand side of Eq. 5.13a. 


definition of an expectation and with the use of Eqs. 5.15 and 5.17, the 
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expected value may be written 


ii -f « | Vy. #, W) di 
Eyal ag 1a Cu. toy V > f, (v ) dv 


where 
pe 








f, (v) = e a (5.19) 


is the probability density function for a zero mean Gaussian random variable 
with unit variance. Expansion of Eq. 5.18 leads to an expression for the 
expectation in terms of eight integrals. 
=. - 
ne are ee a ® 
E eae a Uo 


~ nj 


=_ OO 
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onj Ving 
nj 
r ue : ; 
+ | oy we f, (v") dy (5.20) 
et nj 
nj 
Let 
nJ 

F, (a, 5) = [ f, (&) dé (5.21) 


where f, (&) is defined by Eq. 5.19. The function F, is simply the proba- 
bility distribution function for a Gaussian random variable with zero mean 
and unit variance. Because of the symmetry of the Gaussian probability 


density function about the mean value, 


-) (5.22) 


Fy (-a nj 


al! = ] =e Fy (a 


Equations 5.21 and 5.22 are helpful in evaluating the integrals on the right 
hand side of Eqs. 5.20. 

The first and fifth integrals in Eq. 5.20 follow directly from 
ease 5.cl and 5.22. 


; . - fy (v) dv = - ae [1 - Fy(on5)J (5.23a) 





eeeeaeeeEeEeEeEeEeEEeEeEeEeEeE—J~~ OE Eee 
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co 


[ uP fy (vy) dv 3 F, (a,.) (5.23b) 
nj nj J 


ii 
fm ) 
Oo 


Ge. 
nJ 


The second and sixth integrals in Eq. 5.20 are determined by direct integra- 


tion and by the use of Eq. 5.19. 


-¢, 3 
ny, ox ok 
- 3 f | uc oF ve fy (v ) dv 
ee oy (3 
= uf eee Cre (5. 24a) 
nj nj J 
- oe 
aa Ma ay v t, (v7) dv 
J nj nj 
nj 
8 ih oy fy (a,.) (5. 24b) 
nj ong J 


By employing the method of integration by parts, the third and seventh in- 


tegrals in Eq. 5.20 are evaluated as 


em 2 «xe sa ae 
- 3 i U oF Vv ty edy 
eae no 


2 
= - 3u O° [aso t: (a .) + | =9F (ow) | (5.25a) 
af Vag ng ll vn ee (0 
r 2 
3 ji i og 0 FCW) 
s nj nj 
nj 
© 2 
os é) U oy LF (a5) ame “ng t (a4) (5.25b) 





Uses 


Finally, the fourth and eighth integrals are given by 


nj «3 ; 
= ae vi - f (v") dv 
V. ] 
ar nJ 
i 3 2 
nJ 
and 
r 3 e ox 
i oF va : f, (V) dv 
: nj 
Sing 
ee Ce ran (5.26b) 
Vag ] nj nj E 


Substitution of Eqs. 5.23 through 5.26 into Eq. 5.20 and rearranging terms 


leads to 
me |] = 20. Cis + ace J ores 
J J nj nj nj J 
°2 2 
+ 2 G6 [u +2o, J] f, (a_.) 
Vad Chg Yj eng 
- a, (We +309 ] (5.27) 
nj nj nj 


The expectation term on the right hand side of Eq. 5.13b remains 
to be evaluated. Because of Eqs. 5.15 and 5.17, this expected value may be 


written as 
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E LV |Vng 


{ 
i} 
=> 
cu 
_ 
ce 
oO 
= 
eas 
+ 
Qa 
= 
cs 
< 
+ 
PO 
a 
S 
-_ 
<e 
—— 
oO 
< 


Expansion of Eq. 5.28 leads to six integrals of the same form as those 
evaluated in Eqs. 5.23 through 5.25. Substitution of these integrals into 
Eq. 5.28 yields 


£ = _ HY, 
a ete aad =) enue 
Cnj Vnj 


2 
+ Me [2 F, (on) - 2 oer f, (a,;) - 1] (5.29) 


When coefficients of F, and fy are combined, there results 


£ me : a2 ve 
E (yglWgll = 2 0G + of 1 Fy loys) 
ea oe fig) 
Chi Vag ] nj 
- [ue +06 ] (5.30) 
nj nj 


Substitution of Eqs. 5.27 and 5.30 into Eqs. 5.13 leads to two 
Simultaneous equations for the two coefficients c and c 5° Solution of 
nj 
these equations results in the following expressions for the optimized co- 


efficients, c and. c.:. 
enj nj 





ey 


“co = Cngsni ~ Sning ~ %o t Sn 5 (a, 3) 
+ 4 Si f, (a5) ~ 2 fee (5. 31a) 
aC, (lof : iG, lle 4 leg) = 
2. 2 ion te, f, (a,5)1 (5.31b) 


It 1s readily shown that Eqs. 5.31 do, in fact, minimize the mean square 


error because the second derivatives of the mean square error are all posi- 








tive. 
3 2 2 
E [ee2]) = 2/6 oa gad (5. 32a) 
aa nj nj 
oC, 
nj 
2 
ee fee.) = 2 (5. 32b) 
ace iy 
nj 
and 
2 ea) 2 
sc _.9 ng nj 
nj “e.. 
nJ 
= 87 Un > 0 (5732c) 
nj 
Tt is interesting to note the behavior of the coefficients Co 
A nJ 
and CH; as the current approaches zero velocity. As Wee OF 
; nj 
U 
ee 
nj O-, 0 
nj 
1 
NS ae 
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and 
] 
f es 
a vV 2 
Therefore, as u. > 0, 
nj 
Ca > Cs - 4 os = C. E oy (5.33a) 
nj nj ven nj 
and 
== 0 (5. 33b) 


Equations 5.33 are identical to the optimized coefficients for the case of 
no current given by Eq. 4.24. 
With the mean square error minimized, the random error vector is 


now dropped from Eq. 5.6 leaving the approximate relationship, 


c. IVI} + [ec] Vy = te} + Co] (5.34) 
Substitution of Eq. 5.34 into Eq. 5.5 yields 
[mJ U3 + fc] +0.) + [Kk] tu 
= Emi + Cc] vi + IK) + (0) (5.35) 


Subtraction of [c] {u.} from both sides of Eq. 5.35 and utilization of the 


relationships 
ae ee ee te 


ieee tt ee ay | 


and 
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n= (are 
transforms Eq. 5.35 into the following. 
[Pied ee Well seeses TL) A 
Sc tN) ec) 10) ce ue) eric) snes 
Equation 5.36 is the linearized version of the equation of motion, Eq. 5.2. 
5.2.3 Calculation of Response Statistics 


By definition, the random riser deflection, Y (x,t), is measured 
from the undeflected position of the riser, which is the x axis. The solu- 
tion of Eq. 5.36 1s simplified somewhat if a transformation of coordinates 


meeused. Let 
Meike) = ¥ Got) = WG) (5.37) 


where by (x}, the mean deflection of the riser, is constant for all time be- 


Cause the response is a stationary random process. It follows that 


* 


Y (x,t) (5. 38a) 


eix.t) 
and 


¥ (x,t) = Y. (x,t) (5. 38b) 


By using Eqs. 5.37 and 5.38, Eq. 5.36 may be rewritten as 
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08% é * 
im] 4Y ified i tee SKI eaters 
= eee [c,~]{U} + [~¢ ~]iui} + {c} (5.39) 


a 


Equation 5.39 contains both time dependent and steady terms. For 
the equation to be valid for all time, the steady terms on each side of 


Eq. 5.39 must be equa!. Thus 
[K] fuy} = ec] tut + {c} (5.40) 


Premultiplication of both sides of Eq. 5.40 by cK]! leads to a solution 


for the mean deflections at the nodes. 
t ae = 7} oss. ee : > 
tuys = [kK] ( (e,~J] tu} + te} ) (5.41) 


Because of Eq. 5.40, the time dependent terms of Eq. 5.39 are re- 
lated by 


ml) + (clita + Ik} iv} 


= oe foe bee) Ay (5.42) 


The solution to Eq. 5.42 may be obtained by analogy to the solution of Eq. 4.31 


which is derived in Chapter 4. From Eq. 4.104, it follows that 


14] 


es [ae BE aay 
mmr tio bBo oe 
ica cates NS ica ce ee 
qa) auae ) Sin (2) (5.43) 


* 
Because Y (x,t) is, by definition, a zero mean random variable, the vari- 


ance is given by 


8 


2 
= S 2) da 5.44 
oye i wy (2) (5.44) 
ni 0 ni ni 


The variance of Nee (tjeis 


Se) sm, 7 


5° 
Y 2 a 
ni ni 


. =e 
= EY] 
2 
= Oy (5.45) 


ni 


From Eqs. 5.44 and 5.45, it follows that 


< PM 


ni 


: [ Sx ve (0) do (5.46) 
0 


Vai ni 


Equations 5.41 and 5.46 may be used to calculate the mean and variance of 


the response. Because the input is assumed to be Gaussian and the problem 
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has been linearized, the response is also Gaussian. Therefore, the probabil- 
ity density function of the response is uniquely determined by the response 


mean and variance, and the solution is complete. 
5.3 Riser Response to Random Top Offset 


Even in the absence of external forces, a marine riser will undergo 
a deflection caused by its own weight whenever the top is offset from a point 
vertically above the bottom. Unless the surface support platform is fixed, 
there is almost always some top offset produced by the lateral motion of 
the platform in response to the exciting forces of wind, wave, and current 
and the restoring forces of a mooring system or a dynamic positioning system. 
While a complete study of the motion of surface support platforms is beyond 
the scope of this thesis, an elementary model for top offset may be used to 
calculate the riser response to the combination of wave and current forces 
and top offset. It will be assumed that the top offset is a random process, 
Bop)» which changes so slowly that the associated riser response is essen- 
tially static. This assumption limits the study to platforms and sea condi- 
tions for which the wave-induced platform motion has a negligible effect on 
the riser behavior. 

The random static deflection produced by a random top offset may 


be expressed simply as 


tY (t)} = tb} ¥, 


aw 


ap (5.47) 


where elements of tb}, the vector of influence coefficients, may be deter- 


mined from Eqs. 2.11 and 2.13 as 


> LK, 


* * -]| 
Bn L ni,n-1 anes Th Wi) K ] (5.48) 


ni,n 
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The random bottom rotation Q, (t) caused by top offset is given by 


nese 2) Shel 
gn (t) = BET y(t) = Bet by yg Ct) (5.49) 








and the random bending moments {M (t)} caused by top offset are determined 


after the manner of Eq. 3.96. 
iM (t)} = Menor ae): (5250) 


Equation 5.47 may be used to express the correlation function of 
the deflections at two nodes in terms of the autocorrelation function of 


the top offset. 


R (ey =) (De bark (zr) (525°) 


ee | ni “nj Y 


ani-nj 


A Fourier transformation of both sides of Eq. 5.51 leads to 


By (2) = by bh, Sy (a) (5.52) 


{ni ing 
If the top offset has a zero mean value, then integration of Eq. 5.52 gives 
the variance of the node deflections in terms of the variance of the top 


offset. 


04 peas (5.53) 
ani top 


In a similar manner, it can be shown that 


by 6 (5.54) 


a 22 
%0 L Ytop 


and 
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n n 


2 2 
O = ae ea oO (5.55) 
Mn “es aN ni mk nk one Y top 


If the top offset has a nonzero mean value, uy , then the mean 
top 
responses may be calculated by adapting Eqs. 5.47, 5.49 and 5.50. 





{uot = 4D} up (5.56) 
: i Yop 
n+] 
u = b, u (5.57) 
a L ~ ‘top 
ie Seed] buene (5.58) 
" Y top 


Figures 5.1 and 5.2, which were derived from the results of Fischer 
and Ludwig,’ summarize the behavior of the mean and variance of the bottom 
rotation caused by top offset. In order to maintain a given bottom rotation, 
the ratio of top offset to riser length must decrease as length ratio in- 


creases. 
5.4 Total Response to Current, Waves, and Top Offset 


Let Y (x,t) be the total response of the riser to current, waves, 


and top offset. 
TUES) a Giga) Gantt) (5.59) 


The response to current and waves has been shown to be Gaussian with mean 
BY and variance oy given by Eqs. 5.41 and 5.46. An elementary model for 
the total riser response may be formulated by assuming that the top offset 
is Gaussian (not necessarily zero mean) and that the random response caused 


by waves and current, Y (x,t), and the random response caused by top offset, 
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(x,t), are independent. A Gaussian top offset implies a Gaussian response, 


>< 


(x,t), because the system is linear as shown by Eq. 5.47. If Y (x,t) and 


(x,t) are Gaussian and independent, then Y (x,t) is also Gaussian?’ with 


>< 


mean value 


u + (5.60) 
ni Yn Yai 


ma 
and variance 
o- = 04 +o (5.61) 
In this elementary model, the effect of top offset is completely described 


by the mean and variance of the top offset. 
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Chapter 6 


SAMPLE PROBLEMS 


6.1 General 


The theory derived in Chapters 4 and 5 was used to study the be- 
havior of a typical riser system under various random environmental condi- 
tions. To make this study, two computer programs were developed and coded 
in FORTRAN IV. The first program, RISER DYNAMIC III, calculates the response 
of a riser to random waves only, and the second, RISER DYNAMIC IV, calculates 
the response to a combination of random waves and a steady, deterministic 
current. The University of Illinois IBM System 360-75 computer with the 
FORTRAN IV, Level G compiler was used to obtain solutions. 

The following input parameters were held constant. 


D = 20 inches 


I = .0621 feet* 

Ee = 29 x 10° pounds per square inch 
w = 248 pounds per foot 

Cy) = 1.0 

Cy 20155 


A finite difference model having 31 equally spaced nodes was used to rep- 
resent the riser. Internal friction damping was taken as two percent of 
critical damping for all calculations. 

Input parameters which were varied include the water depth, top 


tension ratio, wind velocity, sea surface elevation spectral density function, 
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and current velocity profile. The natural modes used in each solution in- 
cluded all modes whose frequencies lie in the range where the input sea sur- 
face elevation spectral density function is significant. 

Both computer programs use an iterative procedure which continues 
until [oy ~J and [055] converge. Convergence was considered to occur when 
the assumed value of oy was within five percent of the calculated value 
for every node and the eed value of 05 5 was within ten percent of the 
calculated value for al! modes peeiaered. . most cases, these convergence 


criteria were setisfied after two to four cycles. 
6.2 Riser Response to Random Waves Alone 


The marine riser response to random waves alone was investigated 
for water depths of 609 feet, 1015 feet, 1523 feet, and 2030 feet and sea 
surface elevation spectra given by the Pierson-Moskowitz!© formula for wind 
velocities of ten to forty knots. For this series of problems a top tension 
ratio of 1.2 was used. Two additional studies were made of the 609 foot 
riser, one with the top tension ratio changed to 1.0, the other with the 
sea surface elevation spectrum given by the Sverdrup-Munk-Bretschneider! > 


formula. 
6.2.1 Bottom Rotation 


Figure 6.1 shows the effect of wind velocity, top tension ratio, 
water depth, and sea surface elevation spectral density formula on the root 
mean square bottom rotation. 


As expected, o. increases with wind velocity. The rate of in- 


0 
O 
crease changes rather abruptly at certain wind velocities and remains fairly 
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constant between these changes. This abrupt change in slope occurs whenever 
the low frequency face of the sea surface elevation spectrum approaches a 
natural frequency of the riser. This behavior is best illustrated in the 
case of the 609 foot riser where pronounced changes in slope occur as the 
fundamenta! mode begins to participate in the response. As shown by Fig. 
6.2, for an input SMB sea surface elevation spectrum, the fundamental mode 
1s beginning to partictpate significantly in the response when the wind 
velocity is thirteen knots. The bottom rotation spectra for the 609 foot 
riser with input PM spectra for twenty, twenty-five, and thirty knot winds 
are shown in Figs. 6.3, 6.4 and 6.5. The participation of the fundamental 
mode in the response is seen to correlate with the increased slope of the 
corresponding curves in Fig. 6.1. Once the fundamental mode has been 
picked up, the slope of the 609 foot riser curves remains fairly constant. 
This slope does not appear in the curves for the longer risers because the 
first mode does not participate significantly in their response for the 
wind velocities studied. 

Figure 6.1 shows that the bottom rotation decreases as riser length 
increases, all other parameters remaining constant. In a mathematical sense, 
this results from the fact that the bottom rotation of the normalized mode 
shapes decreases as riser length increases. In a physical sense, the reason 
for this behavior is that the random wave forces become relatively closer 
to the upper end of the riser as length increases. 

The decrease in bottom rotation with increasing water depth is 
not the result of increased damping. On the contrary, damping ratios de- 
crease somewhat as water depth increases, as shown in Fig. 6.6. This be- 


havior may be understood by recalling that hydraulic damping is a function 
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of the relative fluid structure velocity. As water depth increases, a de- 
crease in Structural response results in smaller relative velocities, as 
shown in Fig. 6.8. 

In Fig. 6.1, the response of the 609 foot riser with top tension 
ratio of |.2 to the 20 knot PM spectrum is somewhat larger than a smooth 
curve through the points for other wind velocities would indicate. The ex- 
planation for this perturbation 1s that the peak frequency of the twenty 
knot PM spectrum and the second mode frequency of the 609 foot riser nearly 
coincide, as shown in Fig. 6.3. Similar perturbations might appear in the 
other curves of Fig. 6.1 1f additional data points were generated where 
natural frequencies of the riser coincide with peak frequencies of the sea 
Surface elevation spectrum. 

The effect of top tension ratio on the bottom rotation of the 
609 foot riser is shown in Fig. 6.1. The reduction in rotation which re- 
sults from increased tension may be traced to two factors. First, increasing 
top tension increases tne natura! frequency somewhat and therefore reduces 
the value of the frequency response function. An increase in top tension 
also decreases the bottom rotation of the normalized mode shapes. 

Comparison of the response of the 609 foot riser to SMB and PM 
Spectra iilustrates the significant effect which the input spectrum has on 
the response. Of course, part of this difference can be attributed to the 
fact that, for a given wind velocity, the SMB spectrum implies a larger 
energy density than the PM spectrum, as shown by Fig. C.4. To eliminate this 
source of difference, So 
though the difference eee the two response curves is smaller than in 


1s plotted versus energy density in Fig. 6.9. Al- 


Fig. 6.1, it 1§ still significant, particularly for wind velocities of 25 
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to 30 knots. Thus the response of the riser 1s dependent not only upon the 
energy dens’?ty of the sea but also upon the distribution of that energy with 
respect to frequency. A good :!lustration of this distribution difference 
1S given in Fig. 6.10, which shows the SMB and PM spectra for an energy den- 
sity of 1000 foot pounds per square foot. Although the areas under these 
spectra are equal, the SMB spectrum is concentrated near the fundamental 
riser frequency, and the PM spectrum is spread over higher frequencies. 

The variation of bottom rotation with energy density shows an 
interesting trend when plotted with energy density on a natural scale as in 
Fig. 6.11. The decrease in initial slope with increasing energy density is 
a result of increased hydraulic damping, as shown in Fig. 6.7. Above an 
energy density of 1000 foot pounds per square foot, the variation of Oy 


O 
= * Si 
with E 1s near'y linear. 


6.2.2 Bending Moments 


The variation of root mean square maximum bending moment with 
wind velocity 1s shown in Fig. 6.12, A comparison of Figs. 6.12 and 6.1 
Shows that the maximum bending moment behaves very much !1ke the bottom ro- 
tation as the prob'iem parameters are varied. However, top tension 1s some- 
what more effective in reducing bending moments than it is in reducing bot- 
tom rotations. While the bending moments induced by the SMB spectrum are 
greater than those caused by the corresponding PM spectrum, the difference 
is not as great as the difference between bottom rotations produced by these 
two spectra. 

For the riser cross section studied, the relationship between bend- 


Ing moment, M, and flexura’ stress, f, is 
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f (kips/square inch) = M (OGHe s) (6.1) 


6.2.3 Number of Modes in Response 


It has already been shown, in Figs. 6.2 through 6.5, that two or 
three natural modes may participate in the response of the 609 foot riser. 
As riser length increases, the natural frequencies decrease and become more 
closely grouped with respect to frequency, as shown in Figs. 4.3 and 4.4. 
Therefore, as water depth increases, more modes and higher modes participate 
in the response. This behavior is illustrated in Figs. 6.13 and 6.14 for 
the 2030 foot riser. The response to the twenty knot PM spectrum contains 
modes three through six, and the second through fifth modes contribute sig- 
nificantly to the response to the thirty knot PM spectrum. Figures 6.15 
and 6.16 show the modal contributions to profiles of root mean square de- 


flection and moment produced by the thirty knot PM spectrum. 
6.3 Riser Response to Random Waves and Steady Deterministic Current 


To investigate the effect of steady current on riser response to 
random waves, a series of studies was made using the 609 foot riser with a 
top tension ratio of !.2 and the 30 knot PM wave spectrum. Two current 
profiles were considered, a uniform current and an exponentially decaying 
Current given by 

= ar (hex) 

. a f 

Un (x) = Un (h) © e (G52) 


where De 1s the depth of frictional influence. Equation 6.2, 1s Ekman's 


equation for the magnitude of a pure wind driven current.?° Although Ekman's 


[ee 


solution indicates that the current direction changes with depth (the Ekman 
spiral), for the present investigations the current was considered to be 


unidirectional. The depth of frictional influence was taken as 200 feet. 
6.3.1 Uniform Current 


Figure 6.17 summarizes the behavior of the bottom rotation as 
uniform current velocity 1§ varied. As expected, both the mean response 
and static response to current alone increase as current velocity increases. 
The static response becomes a better predictor of the mean response as cur- 
rent velocity increases. 

As current velocity increases, the standard deviation of the bot- 
tom rotation decreases somewhat because of rather large increases in the 
damping ratios, as shown in Fig. 6.!8. These increased damping ratios may 
be traced to increased hydraulic damping coefficients at the nodes of the 
riser. Fig. 6.19 shows the ratio of hydraulic damping coefficients for a 
uniform current of one-half knot to the same coefficients for no current. 
The relative velocity profile of Fig. 6.20 helps explain the behavior of 
the hydraulic damping coefficients. The uniform current causes the great- 
est relative change in hydraulic damping near the midpoint and lower end 
of the riser where the relative velocity 1s sma!'est. Where the relative 
velocity 1s large with respect to the steady current, the current has little 
effect on hydraulic damping. 

Figures 6.2! and 6.22 show the effect of a one-half knot uniform 
Current on defiection and bending moment statistics of the 609 foot riser 
with a thirty knot PM spectrum. The mean deflections and moments are some- 


what greater than the static response to current alone, and near the top of 
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the riser the mean moment is considerably greater than the static moment. 
Increased hydraulic damping results in standard deviations of deflection 


and moment which are lower with the current than without it. 
6.3.2 Wind Driven Current 


The behavior of the bottom rotation as wind driven current is 
varied is summarized in Fig. 6.23. The increase in mean and static rotations 
is small with respect to the standard deviation, and on a relative basis 
the mean rotation is considerably greater than the static rotation. Unlike 
the uniform current, the wind driven current produces a slight increase in 
the standard deviation of bottom rotation and only a small increase in damp- 
ing ratios for the first two modes, as shown in Fig. 6.24. The reason for 
this behavior is evident from the velocity profiles of Fig. 6.25. The wind 
driven current increases the water velocity in the upper half of the riser 
where it is already large as a result of waves and has a negligible effect 
in the lower half. Therefore the hydraulic damping coefficients are increased 
Slightly in the upper half of the riser and unaffected in the lower half, 
as shown by Fig. 6.26. 

Figures 6.27 and 6.28 show the effect of a wind driven current 
with surface velocity of two knots on the deflection and moment statistics. 
The current causes a slight increase in the standard deviations, and the 


mean response 1s nearly twice the static response. 
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Chapter 7 
CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER STUDY 
7.1 Conclusions 


In this thesis a mathematical model has been developed for pre- 
dicting the nondeterministic response of a marine riser to random wave 
forces, current forces, and random top offset. In this model the continu- 
ous riser structure is represented as a finite degree of freedom system by 
means of the method of finite differences. The validity of this structural 
model and the lumped-smeared force model used to represent wave force dis- 
tributions is demonstrated in Chapter 2. In Chapter 3 a method is derived 
to calculate the static riser response to random waves, and in Chapter 4 
the method is extended to the determination of the dynamic riser response 
to random waves. Chapter 5 completes the derivation of the mathematical 
model by introducing procedures for calculating the response to a combina- 
tion of random waves, a Steady deterministic current, and random top offset. 
In Chapter 6 the mode! is used to demonstrate the effect of several problem 
parameters on riser response. 

From the derivations and sample problems of the preceding five 
chapters several conclusions may be drawn. 

1. The lumped-smeared force model accurately represents wave 
force distributions with a relatively small number of nodes, 
even when the wave force is essentially concentrated at the 
upper end of the riser. 

2. The finite difference structural model, when used with the 


lumped-smeared force model, is a satisfactory representation 
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of the marine riser for purposes of analyzing the riser's be- 
havior in its operating environment. 

The mathematical model derived herein permits an estimation 
of the random dynamic structural response of marine risers 
caused by random waves, steady currents, and random top off- 
sets. 

The high frequency tail of the input sea surface elevation 
Spectrum contributes little to random water velocity and ac- 
celeration. 

The dynamic response of marine risers to random wave forces 
1s significantly greater than the static response. 

The response to random wave forces is highly dependent upon 
the spectral distribution of the wave energy density as well 
as the total average energy density. 

The dynamic response of the riser to waves which are described 
by a given sea surface elevation spectrum decreases as water 
depth increases. 

The random dynamic bottom rotation caused by waves is reduced 
somewhat by increasing the top tension. 

The dynamic response of the riser to random waves may include 
Significant participation by more than one natural mode, the 
lowest one of which is not necessarily the fundamental mode. 
The number of modes which contribute appreciably to the re- 
sponse increases as water depth increases. 

Wave structure interaction results in modified drag forces 
and hydraulic damping which increases significantly with in- 


creased energy density of the sea. 
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For a current velocity profile which is small with respect 
to the standard deviation of the relative water velocity, 
the effect of current on the riser bottom rotation is smal] 
and may be neglected without appreciable error. 

For a current velocity profile which is large with respect 
to the standard deviation of the relative water velocity, the 
current increases hydraulic damping significantly and thus 
reduces the response variance. When this condition exists, 
the static response to current is a good approximation of 
the mean response to current and waves, and the response 
variance for waves alone is a conservative estimate of the 


response variance for the combination of waves and current. 


7.2 Recommendations for Further Study 


The model derived herein 1s somewhat limited by the assumptions 


which were made in its derivation as wel] as the scope of this thesis. 


Several areas where further research would improve the model may be cited. 


ile 


The modei should be tested by comparing calculated riser be- 
havior with field observations of actual riser behavior in 
the ocean environment. A monitoring device for the bottom 
rotation is presently operational*® and could provide data 
which, along with synoptic wave, current, and top offset ob- 
servations, could be used to verify the model. 

Further studies of the platform motion are recommended. At- 
tention should be given to the effect of wave induced plat- 


form motion on riser response. Spectral density functions 
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and probability distribution functions for platform motion 
need to be determined to improve the model. The dynamic ef- 
fect of operational platform motion in deep water, where the 
fundamental frequency may approach platform motion frequencies, 
should also be investigated. 

Because of the cyclic nature of the loads imposed on the riser 
by its environment, fatigue failure is of concern. It is rec- 
ommended that the model developed herein be extended to the 
study of fatigue in marine risers. 

Better information is needed on the spectral distribution of 
the energy of a random sea. While existing wave spectra may 
be satisfactory for wave forecasting, it has been demonstrated 
herein that different spectra lead to greatly different re- 
sponses because the riser is very sensitive to the spectral 
distribution of wave energy. As good directional wave spectra 
are developed, the riser model should also be extended to 
three dimensions. 

The drag and inertia coefficients used in the force formula 
are based on experiments with relatively rigid piles. More 
information is needed on drag and inertia coefficients which 
are appropriate for the rather flexible marine riser. 

The method of analysis employed in this thesis requires the 
use of a digital computer. A simple method should be developed 
for estimating the dynamic response of the riser using hand 


computations only. 
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Tablerceal 


COMPARISON OF MARINE RISER NATURAL FREQUENCIES 


ry ho 3 
ag Finite Finite Finite 
G 7 Series Difference Series Difference Series Difference 
Solution Solution Solution Solution Solution. Solution 
0 3.1416 351403 6.2832 6.2731 9.4248 9.3908 
200 5.4652 5.4638 8.4625 8.4526 iles084 e273 
400 6.3076 6.3061 9.5456 9.5347 12.4565 ize23s 
600 6.8804 6.8789 iOesis 10.2994 13.3145 13.2796 
800 7 3261 7.3245 10.9171 10.9044 14.0089 13.9722 
1000 7.6957 7.6942 11.4243 112-4110 14.5980 14.5597 
7 
=, 
Gy a1 0 
a 
ee (Et 
L m 


Finite difference solution uses 31 equally spaced nodes. 
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Table 2.2 


COMPARISON OF SERIES AND FINITE DIFFERENCE SOLUTIONS 


r 
— 
me 1.0 
Y 

19 lb/ft 

uniform 

load 

7PYNYTKN 
Method Number Maximum Top Bottom 
of of Moment Rotation Rotation 

Solution Nodes (ft-k) (radians) (radians) 
Series -- 165.7 - .0498 1073 
fhite 5 156.6 - .0484 .09505 
a 8 153.6 - .0493 1030 
equal 
fode 1] eles - .0496 - 1058 
spacing) 17 162.0 - .0499 . 1078 
Finite 
Difference 1] W3e71 -.05] 1134 
(unequal 
node spacing) 17 V677 - .0507 sie 
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Table 2.2 Continued 


i 
ieee 
ala WZ 
Vv 
19 lb/ft 
uniform 
load 
2 Y, 

Method Number Maximum Top Bottom 
of of Moment Rotation Rotation 
Solution Nodes (ft-k) (radians) (radians) 

Series -- 94.7 -.0375 . 0682 
Finite 5 88.9 - .0357 0595 

Difference 

(equal 1] 92.3 - .0370 .0660 
node 

Spacing) hs 93.0 - .0372 .0673 
Finite 1] 96.7 - .0375 .0690 
Difference 

(unequal 17 94.9 - .0376 . 0689 


node spacing) 


igi Sia 


Solution 


Serves 


Finite 
Difference 
(equal 
node 
spacing) 


Finite 
Difference 
(unequal 

node spacing) 





Number 
of 
Nodes 


1] 


1] 
We 


Table 2.2 Continued 


Maximum 
Moment 
(ft-k) 


49.6 


46.4 
48.2 
48.5 


50.6 
49.8 


Top 
Rotation 
(radians) 


.0204 


.0214 
.0207 
.0206 


. 0204 
0204 
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Bottom 
Rotation 
(radians) 


.0756 


.0710 
.0745 
BOE | 


.0761 
07654 
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Table 2.2 Continued 


T 


18 AG 


- 


Soe lie 







linear 
current 
variation 


4.5 1b/ft 

Method Number Maximum Top Bottom 

of of Moment Rotation Rotation 
Solution Nodes (ft-k) (radians) (radians) 
Series -- 3/ 35 - .0229 S032 
Finite 5 S50! -.0211 .0276 
Difference 
(equal 1] 36.4 -.0223 . 0300 
node 
spacing) iW 36.8 -.0226 -0305 
Finite 1] 38.4 -.0227 20212 
Difference 
(unequal 17 379 - .0228 .0311 


node spacing) 
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Table 2.2 Continued 


13 |b/tt 







Lopsotrset = 12 ft 


linear 
current 
variation 


Method Number Maximum Top Bottom 
of of Moment Rotation Rotation 

Solution Nodes (ft-k) (radians) (radians) 

Series -- Yana) -.0127 . 069 

Finite 5 oes -.0104 063) 

Difference 

(equal 1] 60.5 -.0120 .0673 

node 

spacing) iz, 60.9 -.0123 .0681 

Finite ul 63.7 -.0125 0692 

Difference 

(unequal 7 62.6 -.0126 .0692 


node spacing) 
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Fig. 1.1 Typical Marine Riser Installation 
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Fig. 1.2 Free Body Diagram of Deflected Riser 
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Fig. 1.3 Forces on Element of Riser Length 
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Fig. 2.1 Node Numbering Scheme 
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Fig. 2.2 Definition Sketch of Constant Tension Beam 
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Fig. 2.4 Spatial Variation of Inertia Force 











Fig. 2.5 Spatial Variation of Drag Force 
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Fig. 2.6 Definition Sketch of Simple 
Beam Subjected to Inertia Force 


174 


weag a[dWig wo} SOLQeY UOLZDeL Jaq UedspLW 7°2 “BL4 


Uy 























(Ux) YSOd _ 4 


— 





(x4) Ysoo 





























(2/) AV corgey uoiqoatseq uedspiy 


Es 


qix) 


tf 7, | Oey th ath, 


OTe ee 


Fig. 2.8 Schematic Representation of Lumped-Smeared Force Model 
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Fig. 2.9 Sketches Used in Derivation of Lumped-Smeared Force Mode 1 
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Fig. 2.10a End Rotation Ratios, Inertia Force, kh = .01 
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Fig. 2.10b End Rotation Ratios, Inertia Force, kh = 10 
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Fig. 2.0 


End Rotation Ratios, Inertia Force, kh = 100 
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Fig. 2.10d End Rotation Ratios, Inertia Force, kh = 10° 
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Fig. 2.1la Midspan Deflection Ratios, Inertia Force, kh = .01 
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Fig. 2.11b Midspan Deflection Ratios, Inertia Force, kh = 10 
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Fig. 2.12b End Rotation Ratios, Drag Force, kh = 10 
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End Rotation Ratios, Drag Force, kh = 100 
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Fig. 2.13b  Midspan Deflection Ratios, Draq Force, kh = 10 
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Fig. 2.13c Midspan Deflection Ratios, Drag Force, kh = 100 
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Fig. 2.13d Midspan Deflection Ratios, Drag Force, kh = 10 
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Fig. 2.14 Natural Frequencies of vertical Marine Riser, & = 1.0 
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Fig. 2.15 Perturbation Parameter 
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Fig. 3.1 Definition Sketch of Sea Surface Elevation 
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Fig. 3.2 Typical Sea Surface Elevation snectru” 
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Fig. 3.3 water Velocity Soectra 
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Appendix A 


ANALYTICAL SOLUTION FOR THE DEFLECTED SHAPE OF A SIMPLE BEAM 
CAUSED BY A WAVE INERTIA FORCE DISTRIBUTION 
Consider a simple beam subjected to a force distribution which 
is characteristic of the inertia force caused by water waves. The govern- 
ing differential equation of the problem, which is depicted in Fig. 2.6, 
is 


d'y _ cosh (kx) 
FT ~ cosh (kh) (A.1) 


The integration of Eq. A.1 four times results in the following expression 
for the deflected shape of the beam 


5) 


y(x) = A, [cosh (kx) +€, “-+¢ Cat C 


gt Co at Cgx + Cy) (A.2) 


where 
A 
= ! 


(kh) EI cosh (kh) 


and Cys Co; C3, and Cy are constants whose values are determined by the 
boundary conditions of the problem. 

For a simply supported beam, the deflection and the curvature 
at the ends of the beam vanish. Application of the two boundary conditions 


at x = 0 yields expressions for the constants Cy and C,. 
Cy = -] (A.3) 


tee 
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The condition of zero curvature at x = h leads to an expression for C,- 
2 2 ] 
C, = (kh)* [1 - cosh (kh)] (A.5) 
] he 


Finally, for zero deflection at x = h, the remaining constant must be 


(kh)? 
3 





2 
- (1 - Ak) cosh (kh) ] r (A.6) 


Substitution of Eqs. A.3 through A.6 into Eq. A.2 results in 
4 


2 5 
y(x) = ee A 7 cosh (kh) 1) 
EI (kh)”~ cosh (kh) 





2 2 2 2 
pin a(x) cy + Ey - MEM) cosh (kn) (4) 


+ cosh (kx) - 1} (Rae) 


As the parameter kh approaches zero, the load approaches a uni- 
formly distributed load of unit magnitude. By employing the infinite series 
representation of the hyperbolic cosine function, it can be shown that the 


limiting form of Eq. A.7 as kh approaches zero is 


é 


lim y(x) = 58 Ee [x’ - 2x°h + xh?] (A.8) 


kh>0 


Equation A.8 is exactly the solution for the deflection of a simply sup- 
ported beam under a unit uniformly distributed load. 
For large values of kh, a simplified form of Eq. A.7 may be written 


by noting that as kh becomes large, 


SSSR TERY approaches zero, 


(kh)° 
cosh (kh) 


and cosh (kx) 
cosh (kh) 


The resulting simplified expression is 


y(x) 


4 
h 
Fi 


large kh 


] 


2 


[---- = 
6 (kh) 


approaches e 


(kh 


ok (x-h) 


4 


] 


ee 


ay 


(kh 


approaches zero, 


4 


k(x-h) 


e 


X 


h 


] 


)} 


6(kh)* 


(*) 


24] 


X 


(A.9) 


| 
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Appendix B 


ANALYTICAL SOLUTIONS FOR THE DEFLECTED SHAPES OF A CONSTANT 
TENSION BEAM CAUSED BY WAVE DRAG AND INERTIA FORCE DISTRIBUTIONS 


B.| Inertia Force Deflections 


Consider the constant tension beam shown in Fig. 2.2 with a force 
distribution which is characteristic of the inertia force caused by water 


waves. The governing differential equation of the problem is 


cp ty. 7 ay cosh (kx) (B.1) 
Ae dx? cosh (kh 


A general solution to Eq. B.1 is 


y = A, {C, + Cyx + C, cosh (ux) + Cy sinh (ux) 
+ A, cosh (kx)} (B.2) 
where 
2S lL. 
4 EI 
Z n* 
* i Eee | 
Elvecosh (kh [kh] 
ee 
2 A,-| 
Z 
pe 
ey 


and C Co; C35 and Cy are coefficients which are determined by the bound- 


ary conditions of the problem. 
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The four boundary conditions are that the deflection and moment 
(or curvature) at each end of the beam are zero. If the curvature at x = 0 


is zero, then 


C2 a -ADA3 (B.3) 


For the curvature to vanish at x = h, it is necessary that 
ae 
Cy = sinh (uh) [cosh (uh) - cosh (kh) ] (B.4) 
The condition of zero deflection at x = 0 implies that 
(B.5) 


Finally, for the deflection to vanish at x = h, the remaining constant 


must be 


ae 
Cy = sgl - cosh (kh) ] (B.6) 


Substitution of Eqs. B.3 through B.6 into Eq. B.2 yields 


4 2 
2 ~~ pk {1 = {il = cosh (kh)] > 
EI (kh)" cosh (kh) © 
2 KZ se Silly ee 
= 3 > cosh ip ate ——7> [cosh (uh) - cosh (kh)] Sink eh 
ko - Ko 4 
Zc 
1 aes CS kx} (B.7) 
-u 


A simplified version of Eq. B.7 may be written for large values 


h h 


of kh where ek approaches zero and cosh (kh) approaches 1/2e* . The re- 


sulting expression for deflection is 
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; he 3a) seen kin pen 
large kh EI (kh)* h 2 sinh (yh) 
A 
Cee Gn os & _ pek(xh) 4, (B.8) 
3 


It is also of interest to examine the left end rotation of the 
constant tension beam. This rotation is analogous to the bottom rotation 
of a marine riser. The general solution for the left end rotation is ob- 


tained by evaluating the first derivative of Eq. B./7 with respect to x at 


' AJA3 Aguh 
ye(o) = Shee {eoshackh) = 1 + sinh (uh) [cosh (uh) 


cosh (kh) ]} (B.9) 


A similar expression for the left end rotation when kh is large 
1s obtained from the first derivative of Eq. B.8 with respect to x, 


evaluated at x = 0. 


“(o) = paged {] + sh fe(uck)h - 14} (B.10) 
een EI rkhy? sinh (uh) 


B.2 Drag Force Deflections 


The governing differential equation of a constant tension beam 
having an applied force distribution which is characteristic of the water 


wave drag force is 


4 @ 2 
d'y .d°y _ ,-cosh (kx) 
eI at U dx? Leosh (kh) (B.11) 
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The general solution to Eq. B.11 is given by 


y = A, {C, + CaXx + C, cosh (ux) + Cy sinh (px) 
2 7X 2 
1: Ao cosh (2kx) - 2A. (kh) (7) } (Baliz) 
where 
_ ht 1 
oi Malas) ae omnes 7. Sameer rar 
8 (kh)" cosh” (kh) 
- ] 
7° ae 
ee Ns 
and 
2 
_ (k 
Az = () 


In order to satisfy the conditions of zero bending moment and zero 


deflection at x = 0, the coefficients C. and C, must be 


o 
! 


ak 4A,[A, - Ao] (B.13) 
and 


2 


3 (B.14) 


‘i A,L4A, - 1] - 4A 
Because the bending moment at x = h is zero, the coefficient Cy 1s given by 


AA 
= Soa Le - A, cosh (2kh) + (Ay - Ag)cosh (uh)] (B.15) 


Finally, for the deflection to vanish at x = h, the remaining coefficient 


must be 
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a, = E tAoL 4A, - 1][cosh (2kh) - 1] + 2A, (kh)? 


y } (B.16) 


The substitution of Eqs. B.13 through B.16 into Eq. B.1I2 leads to 


2 


y= AjAg t-2 (kh)® (4) + [cosh (2kh) - 1 + 2(kh)*] (%) 


+ [1 - 4A,] + 4 [A, - Ay] cosh (ux) 
+ SETH [A, = Ao cosh (2kh) as (A, = A3) cosh (uh) ]sinh (ux) 


A, 
+ A cosh (2kx) } (Bele) 
3 


Equation B.17 may be simplified somewhat for large values of kh, where 


e72kh - 9 

eosh 2(kh) = rev 
and 

Gosh (2kh) = sen 


The resulting simplified expression for the deflection is 


4 


e h i x (u-2k)h sinh ux 
3 ET waa tAg F - ESE + 4A,A,] 27 * 
We Le A» ange (B.18) 


where C, 1S given by Eq. B.13. 
The left end rotation for the general case is obtained from the 


first derivative of Eq. B.12 with respect to x, evaluated at x = 0. 
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y (0) = A, {C,, + uCy} (B.19) 


An expression for left end rotation when kh is large is obtained 


from Eq. B.18 
3 
io) = Fog - SRG 
large kh A(kh) 
ic e(u-2k)h , 4A5A3]} (B.20) 
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Appendix C 
SEA SURFACE ELEVATION SPECTRA 


Oceanographic literature contains several formulas for the one- 
dimensional sea surface elevation spectral density function. These formulas 
have been derived by applying time series analysis techniques to ocean wave 
records. Each of the formulas differs somewhat from the others because of 
differences in the raw data analyzed as well as differences in methods of 
analysis. 

Three of the more commonly used spectral density functions are 
those developed by Pierson, Neumann, and James!* (PNJ); Sverdrup, Munk, and 
Bretschneider!*> (SMB); and Pierson and Moskowitz!& (PM), each of which is 
defined somewhat differently. While the area under the PM spectrum is equal 
to the variance of the wave record,:® the integral of the PNJ spectrum is 
equal to twice the variance of the wave record, and the integral of the SMB 
Spectrum is eight times the variance of the wave record.!! In this thesis, 
the sea surface elevation spectral density is defined in the same way as 
the PM spectral density function. Thus, the variance of the wave record 1s 


given by 


nn 


= = [ s (a) do (C.1) 
O 


Therefore, the PNJ and SMB spectral density functions given here are one- 
half and one-eighth, respectively, of the spectral density functions usually 
found in the literature. 

The PNJ spectral density function for a fully developed sea, as 


given by Kinsman?® and modified to satisfy Eq. C.1, is 
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5 = 25.8 0°° exp [-2 (-2)°] (C2 
"MPN W 
where _ is in feet“-seconds, 2 1s the circular frequency in radians per 
PNJ 


second, and am is the wind velocity. For a fully developed sea, the SMB 


spectrum,/° modified to agree with Eq. C.1, is given by 


= 0.959 2 exp [- 229) (C.3) 


nn (v2) 
where Saneng is again in feet“ seconds and the wind velocity is in knots. 
The Pierson-Moskowitz spectrum!® is given by 

op 7 8-4 0? exp [-0.74 (E)*] (C.4) 

PM Ww 

Figures C.1, €.2, and C.3 show the three spectra for fully developed seas 
and wind velocities of 20, 30, and 40 knots. For a given wind velocity, 
the shapes of the three spectra are different. In particular, the peak of 
the SMB spectrum is sharper and is located at lower frequencies than the peaks 
of the other two spectra. 

It can be shown that the energy density, or total average energy 
per unit surface area of the sea, is equal to the unit weight of the water 
multiplied by the variance of the sea surface elevation.2© Thus, the area 
under the spectral density function is a measure of the energy density. Inte- 
gration of Eqs. C.2, C.3, and C.4 leads to expressions for the energy density, 
E in foot-pounds per square foot of sea surface, in terms of the wind velo- 


eity, Vy? in knots. 
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is ; -3 4 
Eom = 2.51 x 10 Vw (C26) 
x =A 
ay = 1239 x 10 fe (C27) 


In Fig. C.4, the variation of energy density with wind velocity is shown for 
the three spectra. Except at wind velocities of 18.25 knots and 32.9 knots, 
where the PNJ curve crosses the PM and SMB curves, each of the three spectral 
density formulas predicts a different energy density for the fully developed 


sea at a given wind velocity. 
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Sea Surface Elevation Spectra, 20 Knot Wind 
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Fig. C.4 Average Energy Densities for Three Wave Spectra 
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